From d5f263639eaff6d2aed7765aaac3bce6203d8fb5 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Mon, 22 Apr 2024 06:04:40 +0000 Subject: [PATCH 01/34] inital lmm_perm --- .gitignore | 1 + Project.toml | 16 +- docs/Project.toml | 2 + docs/literate/tutorial/lmm_clusterdepth.jl | 71 +++++++++ ext/UnfoldStatsMixedModelsPermutationsExt.jl | 147 +++++++++++++++++++ 5 files changed, 236 insertions(+), 1 deletion(-) create mode 100644 docs/literate/tutorial/lmm_clusterdepth.jl create mode 100644 ext/UnfoldStatsMixedModelsPermutationsExt.jl diff --git a/.gitignore b/.gitignore index 91a6912..401934b 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,4 @@ /docs/Manifest.toml /docs/build/ /docs/src/generated +/docs/dev diff --git a/Project.toml b/Project.toml index f8841fc..e1ad360 100644 --- a/Project.toml +++ b/Project.toml @@ -5,11 +5,25 @@ version = "1.0.0-DEV" [deps] BSplineKit = "093aae92-e908-43d7-9660-e50ee39d5a0a" +ClusterDepth = "c8d8bbfa-f476-4995-adff-2987f04015d1" +MixedModelsPermutations = "647c4018-d7ef-4d03-a0cc-8889a722319e" StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" Unfold = "181c99d8-e21b-4ff3-b70b-c233eddec679" +[weakdeps] +ClusterDepth = "c8d8bbfa-f476-4995-adff-2987f04015d1" +MixedModelsPermutations = "647c4018-d7ef-4d03-a0cc-8889a722319e" + +[extensions] +UnfoldStatsMixedModelsPermutationsExt = ["ClusterDepth", "MixedModelsPermutations"] + [compat] -julia = "1" +BSplineKit = "0.16,0.17" +ClusterDepth = "0.2" +MixedModelsPermutations = "0.2" +StatsModels = "0.7" +Unfold = "0.6,0.7" +julia = "1.9,1.10" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/docs/Project.toml b/docs/Project.toml index 2e1ad97..2af77c4 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,12 +1,14 @@ [deps] CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Chain = "8be319e6-bccf-4806-a6f7-6fae938471bc" +ClusterDepth = "c8d8bbfa-f476-4995-adff-2987f04015d1" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589" +MixedModelsPermutations = "647c4018-d7ef-4d03-a0cc-8889a722319e" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Unfold = "181c99d8-e21b-4ff3-b70b-c233eddec679" diff --git a/docs/literate/tutorial/lmm_clusterdepth.jl b/docs/literate/tutorial/lmm_clusterdepth.jl new file mode 100644 index 0000000..3aab3b2 --- /dev/null +++ b/docs/literate/tutorial/lmm_clusterdepth.jl @@ -0,0 +1,71 @@ + +# get some data +using UnfoldSim +using Unfold +using MixedModelsPermutations, ClusterDepth # both necessary to activate correct extension! + +design = MultiSubjectDesign(; + n_subjects = 30, + n_items = 40, + items_between = Dict(:stimtype => ["car", "face"]), +) +#both_within = Dict(:condition=>["scrambled","intact"])) +contrasts = Dict(:stimtype => DummyCoding()) +p1 = MixedModelComponent(; + basis = UnfoldSim.p100(; sfreq = srate), + formula = @formula(dv ~ 1 + (1 | subject) + (1 | item)), + β = [5.0], + σs = Dict(:subject => [0.0], :item => [0.0]), + contrasts = contrasts, +); + +n1 = MixedModelComponent(; + basis = UnfoldSim.n170(; sfreq = srate), + formula = @formula(dv ~ 1 + stimtype + (1 + stimtype | subject) + (1 | item)), + β = [1.0, 4], # n170-basis is negative + σs = Dict(:subject => [2.0, 0.25], :item => [0.25]), + contrasts = contrasts, +); + +p3 = MixedModelComponent(; + basis = UnfoldSim.p300(; sfreq = srate), + formula = @formula(dv ~ 1 + (1 | subject) + (1 + stimtype | item)), + β = [4.0], + σs = Dict(:subject => [1.0], :item => [0.5, 2]), + contrasts = contrasts, +); + + + +data_e, events = UnfoldSim.simulate( + MersenneTwister(23), + design, + [p1, n1, p3], + UniformOnset(srate * 2, 10), + PinkNoise(; noiselevel = 1); + return_epoched = true, +) # 18 +times = range(-0.1, 0.5, length = size(data_e, 1)) + +#events.latency .+= repeat(range(0,length=size(data,2),step=size(data,1)),inner=size(events[events.subject.=="S01",:],1)) +#data_e,times = Unfold.epoch(data[:],events,[-0.1,0.5],srate) +#events,data_e = Unfold.drop_missing_epochs(events,data_e) + + +# # Fit LMM +m = fit( + UnfoldModel, + Dict( + Any => ( + @formula(0 ~ 1 + stimtype + (1 + stimtype | item) + (1 + stimtype | subject)), + times, + ), + ), + events, + data_e, +); + + +# # Cluster Permute :) +coefficient = 2 +pvalue(MersenneTwister(1), m, data_e, coefficient; n_permutations = 100) \ No newline at end of file diff --git a/ext/UnfoldStatsMixedModelsPermutationsExt.jl b/ext/UnfoldStatsMixedModelsPermutationsExt.jl new file mode 100644 index 0000000..d7ef0a3 --- /dev/null +++ b/ext/UnfoldStatsMixedModelsPermutationsExt.jl @@ -0,0 +1,147 @@ +module UnfoldStatsMixedModelsPermutationsExt +using Unfold +using UnfoldStats +#using MixedModels +using MixedModelsPermutations +using ClusterDepth +using Logging + +function UnfoldStat.pvalue( + rng, + model::UnfoldLinearMixedModel, + data, + coefficient; + type = "clusterdepth", + kwargs..., +) + if type != "clusterdepth" + error("other types (e.g. FDR currently not implemented") + elseif type == "clusterdepth" + return lmm_clusterdepth(rng, model, data, coefficient; kwargs...) + end +end + + + + +function lmm_clusterdepth( + rng, + model, + data, + coefficient; + lmm_statistic = :z, + clusterforming_threshold, + kwargs..., +) + permuted = lmm_permutations(rng, model, data, coefficient; kwargs...) + observed = get_lmm_statistic(model, coefficient, lmm_statistic) + return lmm_clusterdepth_pvalues( + rng, + observed, + permuted; + clusterforming_threshold, + kwargs..., + ) +end +function lmm_clusterdepth_pvalues(observed, permuted; clusterforming_threshold) + + # we need global variables here (yes, sorry...), because instead of actually + # letting ClusterDepth do the permutation, we just have to index the already + # permuted data given in the function (`permuted`) + global n_permutation_count + n_permutation_count = 0 + function _fake_permutation_fun(r, data) + global n_permutation_count + n_permutation_count = n_permutation_count + 1 + return permuted[:, n_permutation_count] + end + J_tuple = ClusterDepth.perm_clusterdepths_both( + MersenneTwister(1), + abs.(permuted), + _fake_permutation_fun, + clusterforming_threshold; + statfun = x -> abs.(x), + nₚ = size(permuted, 2), + ) + + pvals = ClusterDepth.pvals(observed, J_tuple, clusterforming_threshold) + +end + +function lmm_permutations( + rng::AbstractRNG, + model, + data, + coefficient::Int; + n_permutations = 500, + lmm_statistic = :z, + time_selection = 1:size(data, 2), +) + permdata = Matrix{Float64}(undef, length(time_selection), n_permutations) + LMMext = Base.get_extension(Unfold, :UnfoldMixedModelsExt) + mm_outer = LMMext.LinearMixedModel_wrapper( + Unfold.formula(model), + data[1, 1, :], + modelmatrix(model), + ) + mm_outer.optsum.maxtime = 0.1 # + + chIx = 1 # for now we only support 1 channel anyway + # + #p = Progress(length(time_selection)) + #Threads.@threads for tIx =1:length(time_selection) + #@showprogress "Processing Timepoints" for + Threads.@threads for tIx = 1:length(time_selection) + + # splice in the correct dataa for residual calculation + mm = deepcopy(mm_outer) + mm.y .= data[chIx, time_selection[tIx], :] + + # set the previous calculated model-fit + updateL!(MixedModels.setθ!(mm, Vector(model.modelfit.θ[time_selection[tIx]]))) + + # get the coefficient + H0 = coef(mm) + # set the one of interest to 0 + H0[coefficient] = 0 + # run the permutation + + perm = undef + Logging.with_logger(NullLogger()) do # remove NLopt warnings + perm = permutation( + deepcopy(rng), # important here is to set the same seed to keep flip all time-points the same + n_permutations, + mm; + β = H0, + progress = false, + #blup_method = MixedModelsPermutations.olsranef, + ) # constant rng to keep autocorr & olsranef for singular models + end + + # extract the test-statistic + + permdata[tIx, :] = get_lmm_statistic(permutations, coefficient, lmm_statistic) + + #next!(p) + end # end for + return permdata +end + +function get_lmm_statistic( + permutations::MixedModelFitCollection, + coefficient::Int, + lmm_statistic, +) + [ + getproperty(m, lmm_statistic) for + m in permutations.coefpvalues if String(m.coefname) == coefnames(mm)[coefficient] + ] + +end +function get_lmm_statistic(model::UnfoldLinearMixedModel, coefficient::Int, lmm_statistic) + return get_lmm_statistic(modelfit(model).collection, coefficient, lmm_statistic) + # r = coeftable(m) + # r = subset(r, :group => (x -> isnothing.(x)), :coefname => (x -> x .!== "(Intercept)")) + # tvals = abs.(r.estimate ./ r.stderror) + # return tvals +end \ No newline at end of file From 152b1f6c99a42c0362bb21a5d529c61c0ec3b9cf Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Mon, 22 Apr 2024 07:34:56 +0000 Subject: [PATCH 02/34] Moved UnfoldMixedModels also to new extension --- Project.toml | 4 ++++ ext/UnfoldStatsMixedModelsExt.jl | 17 +++++++++++++++++ ext/UnfoldStatsMixedModelsPermutationsExt.jl | 1 + src/extract_coefs.jl | 8 -------- 4 files changed, 22 insertions(+), 8 deletions(-) create mode 100644 ext/UnfoldStatsMixedModelsExt.jl diff --git a/Project.toml b/Project.toml index e1ad360..68435be 100644 --- a/Project.toml +++ b/Project.toml @@ -6,20 +6,24 @@ version = "1.0.0-DEV" [deps] BSplineKit = "093aae92-e908-43d7-9660-e50ee39d5a0a" ClusterDepth = "c8d8bbfa-f476-4995-adff-2987f04015d1" +MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316" MixedModelsPermutations = "647c4018-d7ef-4d03-a0cc-8889a722319e" StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" Unfold = "181c99d8-e21b-4ff3-b70b-c233eddec679" [weakdeps] ClusterDepth = "c8d8bbfa-f476-4995-adff-2987f04015d1" +MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316" MixedModelsPermutations = "647c4018-d7ef-4d03-a0cc-8889a722319e" [extensions] +UnfoldStatsMixedModelsExt = ["MixedModels"] UnfoldStatsMixedModelsPermutationsExt = ["ClusterDepth", "MixedModelsPermutations"] [compat] BSplineKit = "0.16,0.17" ClusterDepth = "0.2" +MixedModels = "4" MixedModelsPermutations = "0.2" StatsModels = "0.7" Unfold = "0.6,0.7" diff --git a/ext/UnfoldStatsMixedModelsExt.jl b/ext/UnfoldStatsMixedModelsExt.jl new file mode 100644 index 0000000..61426b9 --- /dev/null +++ b/ext/UnfoldStatsMixedModelsExt.jl @@ -0,0 +1,17 @@ +module UnfoldStatsMixedModelsExt +using Unfold +using UnfoldStats + +lmm_ext = Base.get_extension(Unfold, :UnfoldMixedModelsExt) +# Currently, `extract_coefs` is not implemented for mixed-effects models +UnfoldStats.extract_coefs( + model::Union{ + lmm_ext.UnfoldLinearMixedModel, + lmm_ext.UnfoldLinearMixedModelContinuousTime, + }, + predictor, + basisname, +) = throw( + "The `extract_coefs` function is currently not implemented for mixed-effects models.", +) +end \ No newline at end of file diff --git a/ext/UnfoldStatsMixedModelsPermutationsExt.jl b/ext/UnfoldStatsMixedModelsPermutationsExt.jl index d7ef0a3..00309b5 100644 --- a/ext/UnfoldStatsMixedModelsPermutationsExt.jl +++ b/ext/UnfoldStatsMixedModelsPermutationsExt.jl @@ -144,4 +144,5 @@ function get_lmm_statistic(model::UnfoldLinearMixedModel, coefficient::Int, lmm_ # r = subset(r, :group => (x -> isnothing.(x)), :coefname => (x -> x .!== "(Intercept)")) # tvals = abs.(r.estimate ./ r.stderror) # return tvals +end end \ No newline at end of file diff --git a/src/extract_coefs.jl b/src/extract_coefs.jl index 0ac0d2c..ef4c94a 100644 --- a/src/extract_coefs.jl +++ b/src/extract_coefs.jl @@ -170,14 +170,6 @@ function extract_coefs(model::UnfoldModel, predictor, basisname) return coef_subset#::Array{<:Union{<:Missing,<:Float64},3} end -# Currently, `extract_coefs` is not implemented for mixed-effects models -extract_coefs( - model::Union{UnfoldLinearMixedModel,UnfoldLinearMixedModelContinuousTime}, - predictor, - basisname, -) = throw( - "The `extract_coefs` function is currently not implemented for mixed-effects models.", -) """ extract_coefs(models::Vector{<:UnfoldModel}, predictor, basisname) From 638ba0b0e59861bb3020ae1d128de836f2188539 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Mon, 22 Apr 2024 12:47:00 +0000 Subject: [PATCH 03/34] well, it does something now --- Project.toml | 3 + docs/Project.toml | 1 + docs/literate/tutorial/lmm_clusterdepth.jl | 23 +++++-- ext/UnfoldStatsMixedModelsExt.jl | 6 ++ ext/UnfoldStatsMixedModelsPermutationsExt.jl | 70 ++++++++++++++------ src/UnfoldStats.jl | 2 + 6 files changed, 76 insertions(+), 29 deletions(-) diff --git a/Project.toml b/Project.toml index 68435be..b90b012 100644 --- a/Project.toml +++ b/Project.toml @@ -6,8 +6,11 @@ version = "1.0.0-DEV" [deps] BSplineKit = "093aae92-e908-43d7-9660-e50ee39d5a0a" ClusterDepth = "c8d8bbfa-f476-4995-adff-2987f04015d1" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316" MixedModelsPermutations = "647c4018-d7ef-4d03-a0cc-8889a722319e" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +StatsAPI = "82ae8749-77ed-4fe6-ae5f-f523153014b0" StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" Unfold = "181c99d8-e21b-4ff3-b70b-c233eddec679" diff --git a/docs/Project.toml b/docs/Project.toml index 2af77c4..9c18e19 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -8,6 +8,7 @@ Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589" +MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316" MixedModelsPermutations = "647c4018-d7ef-4d03-a0cc-8889a722319e" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" diff --git a/docs/literate/tutorial/lmm_clusterdepth.jl b/docs/literate/tutorial/lmm_clusterdepth.jl index 3aab3b2..8cd5732 100644 --- a/docs/literate/tutorial/lmm_clusterdepth.jl +++ b/docs/literate/tutorial/lmm_clusterdepth.jl @@ -3,7 +3,11 @@ using UnfoldSim using Unfold using MixedModelsPermutations, ClusterDepth # both necessary to activate correct extension! +using UnfoldStats +using StatsModels +using Random +srate = 25 design = MultiSubjectDesign(; n_subjects = 30, n_items = 40, @@ -38,7 +42,6 @@ p3 = MixedModelComponent(; data_e, events = UnfoldSim.simulate( - MersenneTwister(23), design, [p1, n1, p3], UniformOnset(srate * 2, 10), @@ -46,21 +49,20 @@ data_e, events = UnfoldSim.simulate( return_epoched = true, ) # 18 times = range(-0.1, 0.5, length = size(data_e, 1)) - +data_e = reshape(data_e, size(data_e, 1), :) #events.latency .+= repeat(range(0,length=size(data,2),step=size(data,1)),inner=size(events[events.subject.=="S01",:],1)) -#data_e,times = Unfold.epoch(data[:],events,[-0.1,0.5],srate) -#events,data_e = Unfold.drop_missing_epochs(events,data_e) + # # Fit LMM m = fit( UnfoldModel, - Dict( + [ Any => ( @formula(0 ~ 1 + stimtype + (1 + stimtype | item) + (1 + stimtype | subject)), times, ), - ), + ], events, data_e, ); @@ -68,4 +70,11 @@ m = fit( # # Cluster Permute :) coefficient = 2 -pvalue(MersenneTwister(1), m, data_e, coefficient; n_permutations = 100) \ No newline at end of file +pvalues( + MersenneTwister(1), + m, + data_e, + coefficient; + n_permutations = 10, + clusterforming_threshold = 1.8, +) \ No newline at end of file diff --git a/ext/UnfoldStatsMixedModelsExt.jl b/ext/UnfoldStatsMixedModelsExt.jl index 61426b9..e3bda90 100644 --- a/ext/UnfoldStatsMixedModelsExt.jl +++ b/ext/UnfoldStatsMixedModelsExt.jl @@ -1,8 +1,14 @@ module UnfoldStatsMixedModelsExt using Unfold using UnfoldStats +using MixedModels +import StatsAPI: pvalue lmm_ext = Base.get_extension(Unfold, :UnfoldMixedModelsExt) + +if isnothing(lmm_ext) + error("Something went wrong with getting the Unfold UnfoldMixedModelsExt extension") +end # Currently, `extract_coefs` is not implemented for mixed-effects models UnfoldStats.extract_coefs( model::Union{ diff --git a/ext/UnfoldStatsMixedModelsPermutationsExt.jl b/ext/UnfoldStatsMixedModelsPermutationsExt.jl index 00309b5..a718fb0 100644 --- a/ext/UnfoldStatsMixedModelsPermutationsExt.jl +++ b/ext/UnfoldStatsMixedModelsPermutationsExt.jl @@ -1,23 +1,36 @@ module UnfoldStatsMixedModelsPermutationsExt using Unfold +import Unfold: pvalues using UnfoldStats #using MixedModels using MixedModelsPermutations using ClusterDepth using Logging +using Random +const MixedModels = MixedModelsPermutations.MixedModels +LMMext = Base.get_extension(Unfold, :UnfoldMixedModelsExt) -function UnfoldStat.pvalue( + +function Unfold.pvalues( rng, - model::UnfoldLinearMixedModel, + model::LMMext.UnfoldLinearMixedModel, data, coefficient; type = "clusterdepth", + clusterforming_threshold, kwargs..., ) if type != "clusterdepth" error("other types (e.g. FDR currently not implemented") elseif type == "clusterdepth" - return lmm_clusterdepth(rng, model, data, coefficient; kwargs...) + return lmm_clusterdepth( + rng, + model, + data, + coefficient; + clusterforming_threshold, + kwargs..., + ) end end @@ -43,7 +56,13 @@ function lmm_clusterdepth( kwargs..., ) end -function lmm_clusterdepth_pvalues(observed, permuted; clusterforming_threshold) +function lmm_clusterdepth_pvalues( + rng, + observed, + permuted; + clusterforming_threshold, + kwargs..., +) # we need global variables here (yes, sorry...), because instead of actually # letting ClusterDepth do the permutation, we just have to index the already @@ -56,7 +75,7 @@ function lmm_clusterdepth_pvalues(observed, permuted; clusterforming_threshold) return permuted[:, n_permutation_count] end J_tuple = ClusterDepth.perm_clusterdepths_both( - MersenneTwister(1), + rng, abs.(permuted), _fake_permutation_fun, clusterforming_threshold; @@ -71,19 +90,18 @@ end function lmm_permutations( rng::AbstractRNG, model, - data, + data::AbstractArray{<:Real,3}, coefficient::Int; n_permutations = 500, lmm_statistic = :z, time_selection = 1:size(data, 2), ) permdata = Matrix{Float64}(undef, length(time_selection), n_permutations) - LMMext = Base.get_extension(Unfold, :UnfoldMixedModelsExt) - mm_outer = LMMext.LinearMixedModel_wrapper( - Unfold.formula(model), - data[1, 1, :], - modelmatrix(model), - ) + + + Xs = LMMext.prepare_modelmatrix(model) + + mm_outer = LMMext.LinearMixedModel_wrapper(Unfold.formulas(model), data[1, 1, :], Xs) mm_outer.optsum.maxtime = 0.1 # chIx = 1 # for now we only support 1 channel anyway @@ -98,7 +116,9 @@ function lmm_permutations( mm.y .= data[chIx, time_selection[tIx], :] # set the previous calculated model-fit - updateL!(MixedModels.setθ!(mm, Vector(model.modelfit.θ[time_selection[tIx]]))) + θ = Vector(modelfit(model).θ[time_selection[tIx]]) + @debug size(θ) + MixedModels.updateL!(MixedModels.setθ!(mm, θ)) # get the coefficient H0 = coef(mm) @@ -106,21 +126,22 @@ function lmm_permutations( H0[coefficient] = 0 # run the permutation - perm = undef + permutations = undef Logging.with_logger(NullLogger()) do # remove NLopt warnings - perm = permutation( + permutations = permutation( deepcopy(rng), # important here is to set the same seed to keep flip all time-points the same n_permutations, mm; β = H0, - progress = false, + hide_progress = true, #blup_method = MixedModelsPermutations.olsranef, ) # constant rng to keep autocorr & olsranef for singular models end # extract the test-statistic - permdata[tIx, :] = get_lmm_statistic(permutations, coefficient, lmm_statistic) + permdata[tIx, :] = + get_lmm_statistic(model, permutations, coefficient, lmm_statistic) #next!(p) end # end for @@ -128,18 +149,23 @@ function lmm_permutations( end function get_lmm_statistic( - permutations::MixedModelFitCollection, + model, + permutations::MixedModelsPermutations.MixedModels.MixedModelFitCollection, coefficient::Int, lmm_statistic, ) [ - getproperty(m, lmm_statistic) for - m in permutations.coefpvalues if String(m.coefname) == coefnames(mm)[coefficient] + getproperty(m, lmm_statistic) for m in permutations.coefpvalues if + String(m.coefname) == Unfold.coefnames(Unfold.formulas(model))[1][coefficient] ] end -function get_lmm_statistic(model::UnfoldLinearMixedModel, coefficient::Int, lmm_statistic) - return get_lmm_statistic(modelfit(model).collection, coefficient, lmm_statistic) +function get_lmm_statistic( + model::LMMext.UnfoldLinearMixedModel, + coefficient::Int, + lmm_statistic, +) + return get_lmm_statistic(model, modelfit(model), coefficient, lmm_statistic) # r = coeftable(m) # r = subset(r, :group => (x -> isnothing.(x)), :coefname => (x -> x .!== "(Intercept)")) # tvals = abs.(r.estimate ./ r.stderror) diff --git a/src/UnfoldStats.jl b/src/UnfoldStats.jl index ae1a2f8..540dcf5 100644 --- a/src/UnfoldStats.jl +++ b/src/UnfoldStats.jl @@ -1,5 +1,6 @@ module UnfoldStats +using BSplineKit: Reexport using Unfold using BSplineKit using StatsModels @@ -8,4 +9,5 @@ include("extract_coefs.jl") # export functions to extract model coefficients export extract_coefs + end From 67fe65355dddb685ac4759190d30c89a3e5b6cbe Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Mon, 22 Apr 2024 12:47:49 +0000 Subject: [PATCH 04/34] added lmm clusterperm --- docs/make.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/make.jl b/docs/make.jl index d5ef8b5..b5c1186 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -35,6 +35,7 @@ makedocs(; "Home" => "index.md", "Tutorials" => ["Two-stage EEG analysis" => "generated/tutorial/two_stage_analysis.md"], + ["LMM EEG + ClusterDepth" => "generated/tutorial/lmm_clusterdepth.md"], ], ) From eaa207d9e86b8373f1d1820cb7d714bb9a259d88 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Mon, 22 Apr 2024 12:56:47 +0000 Subject: [PATCH 05/34] forgot to save this file --- docs/literate/tutorial/lmm_clusterdepth.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/tutorial/lmm_clusterdepth.jl b/docs/literate/tutorial/lmm_clusterdepth.jl index 8cd5732..98d4e85 100644 --- a/docs/literate/tutorial/lmm_clusterdepth.jl +++ b/docs/literate/tutorial/lmm_clusterdepth.jl @@ -49,7 +49,7 @@ data_e, events = UnfoldSim.simulate( return_epoched = true, ) # 18 times = range(-0.1, 0.5, length = size(data_e, 1)) -data_e = reshape(data_e, size(data_e, 1), :) +data_e = reshape(data_e, 1, size(data_e, 1), :) #events.latency .+= repeat(range(0,length=size(data,2),step=size(data,1)),inner=size(events[events.subject.=="S01",:],1)) From 247623ad9eef96e661115028a65b8630f871cc69 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Mon, 22 Apr 2024 12:57:05 +0000 Subject: [PATCH 06/34] forgot to save --- ext/UnfoldStatsMixedModelsExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/UnfoldStatsMixedModelsExt.jl b/ext/UnfoldStatsMixedModelsExt.jl index e3bda90..632373e 100644 --- a/ext/UnfoldStatsMixedModelsExt.jl +++ b/ext/UnfoldStatsMixedModelsExt.jl @@ -2,7 +2,7 @@ module UnfoldStatsMixedModelsExt using Unfold using UnfoldStats using MixedModels -import StatsAPI: pvalue + lmm_ext = Base.get_extension(Unfold, :UnfoldMixedModelsExt) From 409f647abdb799e8b4b0e874f352f7824bbcd5cb Mon Sep 17 00:00:00 2001 From: Benedikt Ehinger Date: Tue, 23 Apr 2024 12:07:20 +0200 Subject: [PATCH 07/34] Update UnfoldStatsMixedModelsPermutationsExt.jl --- ext/UnfoldStatsMixedModelsPermutationsExt.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/UnfoldStatsMixedModelsPermutationsExt.jl b/ext/UnfoldStatsMixedModelsPermutationsExt.jl index a718fb0..c90be37 100644 --- a/ext/UnfoldStatsMixedModelsPermutationsExt.jl +++ b/ext/UnfoldStatsMixedModelsPermutationsExt.jl @@ -133,7 +133,7 @@ function lmm_permutations( n_permutations, mm; β = H0, - hide_progress = true, + progress = false, #blup_method = MixedModelsPermutations.olsranef, ) # constant rng to keep autocorr & olsranef for singular models end @@ -171,4 +171,4 @@ function get_lmm_statistic( # tvals = abs.(r.estimate ./ r.stderror) # return tvals end -end \ No newline at end of file +end From 767e5b14f8b057cf24c2e82d8544640aa994a5fa Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Tue, 30 Apr 2024 10:50:24 +0000 Subject: [PATCH 08/34] added an 'abs' --- ext/UnfoldStatsMixedModelsPermutationsExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/UnfoldStatsMixedModelsPermutationsExt.jl b/ext/UnfoldStatsMixedModelsPermutationsExt.jl index a718fb0..fe7227a 100644 --- a/ext/UnfoldStatsMixedModelsPermutationsExt.jl +++ b/ext/UnfoldStatsMixedModelsPermutationsExt.jl @@ -83,7 +83,7 @@ function lmm_clusterdepth_pvalues( nₚ = size(permuted, 2), ) - pvals = ClusterDepth.pvals(observed, J_tuple, clusterforming_threshold) + pvals = ClusterDepth.pvals(abs.(observed), J_tuple, clusterforming_threshold) end From d747591aaa81a405b3c99b691d005e7480a9cf2e Mon Sep 17 00:00:00 2001 From: Benedikt Ehinger Date: Wed, 22 May 2024 13:53:25 +0200 Subject: [PATCH 09/34] Update extract_coefs.jl --- src/extract_coefs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/extract_coefs.jl b/src/extract_coefs.jl index 0ac0d2c..8d68cab 100644 --- a/src/extract_coefs.jl +++ b/src/extract_coefs.jl @@ -35,7 +35,7 @@ extract_symbol(t::FunctionTerm) = extract_symbol.(t.args) extract_symbol(t::MatrixTerm) = extract_symbol(t.terms) extract_symbol(t::Unfold.TimeExpandedTerm) = repeat(extract_symbol(t.term), inner = length(Unfold.colnames(t.basisfunction))) -extract_symbol(f::FormulaTerm) = vcat(extract_symbol.(f.rhs)...) +extract_symbol(f::FormulaTerm) = extract_symbol(f.rhs) extract_symbol(t::Vector) = vcat(extract_symbol.(t)...) extract_symbol(t::Tuple) = vcat(extract_symbol.(t)...) From 59f07623f42fbe4d69d11ef6bf39104084b57f4c Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Sun, 23 Jun 2024 21:01:31 +0000 Subject: [PATCH 10/34] upgrade to unfold07 --- benchmark/sim_and_fit.jl | 23 ++++++++++++----------- src/extract_coefs.jl | 38 +++++++------------------------------- test/extract_coefs.jl | 35 ++++++++++++++++++----------------- test/setup.jl | 4 ++-- 4 files changed, 39 insertions(+), 61 deletions(-) diff --git a/benchmark/sim_and_fit.jl b/benchmark/sim_and_fit.jl index aec830b..9a7a9db 100644 --- a/benchmark/sim_and_fit.jl +++ b/benchmark/sim_and_fit.jl @@ -36,7 +36,8 @@ function sim_and_fit( ) # At the moment, the function is not implemented for mixed models - if model_type in [UnfoldLinearMixedModel, UnfoldLinearMixedModelContinuousTime] + ext = Base.get_extension(Unfold, :UnfoldMixedModelsExt) + if model_type in [ext.UnfoldLinearMixedModel, ext.UnfoldLinearMixedModelContinuousTime] throw("Not implemented.") end @@ -117,8 +118,8 @@ function create_event_dict( t_stim = times t_fix = times else # UnfoldLinearModelContinuousTime - t_stim = firbasis(τ = (-0.1, 1.5), sfreq = 100, name = "stimulus") - t_fix = firbasis(τ = (-0.1, 1), sfreq = 100, name = "fixation") + t_stim = firbasis(τ = (-0.1, 1.5), sfreq = 100, name = "stim") + t_fix = firbasis(τ = (-0.1, 1), sfreq = 100, name = "fix") end # Define formula(s) @@ -144,8 +145,8 @@ function create_event_dict( t_stim = times t_fix = times else # UnfoldLinearModelContinuousTime - t_stim = firbasis(τ = (-0.1, 1.5), sfreq = 100, name = "stimulus") - t_fix = firbasis(τ = (-0.1, 1), sfreq = 100, name = "fixation") + t_stim = firbasis(τ = (-0.1, 1.5), sfreq = 100, name = "stim") + t_fix = firbasis(τ = (-0.1, 1), sfreq = 100, name = "fix") end # Define formula(s) @@ -153,9 +154,9 @@ function create_event_dict( f_fix = @formula 0 ~ 1 + spl(continuous, 4) # Combine basis function(s)/times and formula(s) with the corresponding event - event_dict = Dict("stim" => (f_stim, t_stim), "fix" => (f_fix, t_fix)) + event_vec = ["stim" => (f_stim, t_stim), "fix" => (f_fix, t_fix)] - return event_dict + return event_vec end function create_event_dict( @@ -176,9 +177,9 @@ function create_event_dict( f = @formula 0 ~ 1 + spl(continuous, 4) + continuous + condition * pet # Combine basis function(s)/times and formula(s) with the corresponding event - event_dict = Dict(Any => (f, t)) + event_vec = [Any => (f, t)] - return event_dict + return event_vec end function create_event_dict( @@ -199,9 +200,9 @@ function create_event_dict( f = @formula 0 ~ 1 + spl(continuous, 4) # Combine basis function(s)/times and formula(s) with the corresponding event - event_dict = Dict(Any => (f, t)) + event_vec = [Any => (f, t)] - return event_dict + return event_vec end get_conditions(sim_type::T) where {T} = get_conditions(PredictorStyle(T)) diff --git a/src/extract_coefs.jl b/src/extract_coefs.jl index ef4c94a..2f16fe5 100644 --- a/src/extract_coefs.jl +++ b/src/extract_coefs.jl @@ -35,7 +35,7 @@ extract_symbol(t::FunctionTerm) = extract_symbol.(t.args) extract_symbol(t::MatrixTerm) = extract_symbol(t.terms) extract_symbol(t::Unfold.TimeExpandedTerm) = repeat(extract_symbol(t.term), inner = length(Unfold.colnames(t.basisfunction))) -extract_symbol(f::FormulaTerm) = vcat(extract_symbol.(f.rhs)...) +extract_symbol(f::FormulaTerm) = extract_symbol(f.rhs) extract_symbol(t::Vector) = vcat(extract_symbol.(t)...) extract_symbol(t::Tuple) = vcat(extract_symbol.(t)...) @@ -63,30 +63,6 @@ julia> UnfoldStats.get_predictor_string(:condition) get_predictor_string(p::Symbol) = ":$p" get_predictor_string(p::String) = "\"$p\"" get_predictor_string(p::Tuple) = "$p" - -""" - get_basisnames(model::UnfoldModel) - -Return the basisnames for all predictor terms as a vector. - -The returned vector contains the name of the event type/basis, repeated by their actual coefficient number (after StatsModels.apply_schema). -If a model has more than one event type (e.g. stimulus and fixation), the vectors are concatenated. -""" -function get_basisnames(model::UnfoldLinearModel) - # Extract the event names from the design - design_keys = keys(Unfold.design(model)) - - # Create a list of the basis names corresponding to each model term - basisnames = String[] - for (ix, event) in enumerate(design_keys) - push!(basisnames, repeat(["event: $(event)"], size(modelmatrix(model)[ix], 2))...) - end - return basisnames -end - -get_basisnames(model::UnfoldLinearModelContinuousTime) = - Unfold.extract_coef_info(Unfold.get_coefnames(model), 1) - """ extract_coefs(model::UnfoldModel, predictor, basisname) @@ -105,7 +81,7 @@ The dimensions of the returned coefficients are channel x times x coefficients. function extract_coefs(model::UnfoldModel, predictor, basisname) # Get vector with underlying predictor variable (symbol) for all coefficients - symbols = extract_symbol(Unfold.formula(model)) + symbols = extract_symbol(Unfold.formulas(model)) # Check whether `predictor` is a predictor in the model if predictor ∉ symbols @@ -119,10 +95,10 @@ function extract_coefs(model::UnfoldModel, predictor, basisname) ) end - basisname_list = get_basisnames(model) + basisname_list = Unfold.get_basis_names(model) # Check whether given `basisname` is in the basisname list of the model - if basisname ∉ basisname_list + if basisname ∉ vcat(basisname_list...) allowed_basisnames = join(["\"$b\"" for b in unique(basisname_list)], ", ") throw( @@ -136,7 +112,7 @@ function extract_coefs(model::UnfoldModel, predictor, basisname) mask_predictor = contained_or_equal.(predictor, symbols) # Create a boolean mask which is true for the model coefficients that belong to the given basis name - mask_basisfunction = basisname_list .== basisname + mask_basisfunction = vcat(basisname_list...) .== basisname mask_combined = mask_predictor .* mask_basisfunction @@ -151,10 +127,10 @@ function extract_coefs(model::UnfoldModel, predictor, basisname) end # Extract the requested coefficients from the coefficient array - if typeof(model) == UnfoldLinearModel + if isa(model, UnfoldLinearModel) coef_subset = coef(model)[:, :, mask_combined] - elseif typeof(model) == UnfoldLinearModelContinuousTime + elseif isa(model, UnfoldLinearModelContinuousTime) n_coefs = length( unique(Unfold.extract_coef_info(Unfold.get_coefnames(model), 2)[mask_combined]), diff --git a/test/extract_coefs.jl b/test/extract_coefs.jl index 6bc20b6..7c855c6 100644 --- a/test/extract_coefs.jl +++ b/test/extract_coefs.jl @@ -8,13 +8,13 @@ σs = Dict(:subject => [1, 0.1, 0.1, 0.1]) simulation = define_simulation(sim_type, β, σs, n_subjects = 5) - + model_type = UnfoldLinearModelContinuousTime for model_type in [UnfoldLinearModel, UnfoldLinearModelContinuousTime] models = sim_and_fit(sim_type, simulation, model_type, seed = 1) model_1 = models[1, :unfoldmodel] - symbols = UnfoldStats.extract_symbol(Unfold.formula(model_1)) + symbols = UnfoldStats.extract_symbol(Unfold.formulas(model_1)) if model_type == UnfoldLinearModel @@ -39,21 +39,21 @@ @test UnfoldStats.get_predictor_string((:condition, :pet)) == "(:condition, :pet)" - @test UnfoldStats.get_basisnames(model_1) == - reduce(vcat, fill.(["event: stim", "event: fix"], [2, 8])) + @test Unfold.get_basis_names(model_1) == + reduce(vcat, fill.(["stim", "fix"], [2, 8])) # Test exceptions - @test_throws ArgumentError extract_coefs(model_1, "continuous", "event: fix") - @test_throws ArgumentError extract_coefs(model_1, :continuous, "fix") - @test_throws ArgumentError extract_coefs(model_1, :pet, "event: stim") + @test_throws ArgumentError extract_coefs(model_1, "continuous", "fix") + @test_throws ArgumentError extract_coefs(model_1, :continuous, "bla fix") + @test_throws ArgumentError extract_coefs(model_1, :pet, "stim") # Test `extract_coefs` method for single Unfold model - coefs_1 = extract_coefs(model_1, :continuous, "event: fix") + coefs_1 = extract_coefs(model_1, :continuous, "fix") @test size(coefs_1) == (3, 15, 4) @test coefs_1[:] == subset( coeftable(model_1), - :basisname => ByRow(==("event: fix")), + :eventname => ByRow(==("fix")), :coefname => ByRow( x -> x in [ "continuous", @@ -65,27 +65,28 @@ ).estimate # Test `extract_coefs` method for a vector of Unfold models - coefs = extract_coefs(models.unfoldmodel, :continuous, "event: stim") + coefs = extract_coefs(models.unfoldmodel, :continuous, "stim") @test coefs[:, :, 1, 4][:] == subset( coeftable(models[4, :unfoldmodel]), - :basisname => ByRow(==("event: stim")), + :eventname => ByRow(==("stim")), :coefname => ByRow(==("continuous")), ).estimate else # Test helper functions @test length(symbols) == 1210 - @test UnfoldStats.get_basisnames(model_1) == - reduce(vcat, fill.(["stimulus", "fixation"], [161 * 2, 111 * 8])) + @test Unfold.get_basis_names(model_1) |> + x -> + vcat(x...) == reduce(vcat, fill.(["stim", "fix"], [161 * 2, 111 * 8])) # Test `extract_coefs` method for single Unfold model - coefs_1 = extract_coefs(model_1, :continuous, "fixation") + coefs_1 = extract_coefs(model_1, :continuous, "fix") @test size(coefs_1) == (3, 111, 4) @test coefs_1[:, 1, :][:] == subset( coeftable(model_1), - :basisname => ByRow(==("fixation")), + :eventname => ByRow(==("fix")), :coefname => ByRow( x -> x in [ "continuous", @@ -98,11 +99,11 @@ ).estimate # Test `extract_coefs` method for a vector of Unfold models - coefs = extract_coefs(models.unfoldmodel, :continuous, "stimulus") + coefs = extract_coefs(models.unfoldmodel, :continuous, "stim") @test coefs[:, :, 1, 4][:] == subset( coeftable(models[4, :unfoldmodel]), - :basisname => ByRow(==("stimulus")), + :eventname => ByRow(==("stim")), :coefname => ByRow(==("continuous")), ).estimate end diff --git a/test/setup.jl b/test/setup.jl index 8a72add..11ca32f 100644 --- a/test/setup.jl +++ b/test/setup.jl @@ -4,5 +4,5 @@ using Unfold using DataFrames using StableRNGs -include("../benchmark/types.jl") -include("../benchmark/sim_and_fit.jl") +includet("../benchmark/types.jl") +includet("../benchmark/sim_and_fit.jl") From 275e4fbf294297fbd8b124c1ab2b73a8f58b831a Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Sun, 23 Jun 2024 21:02:08 +0000 Subject: [PATCH 11/34] remove compatability with unfold 0.6 --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index b90b012..ccc7a70 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "UnfoldStats" uuid = "96fd419a-8306-4ce8-ba5b-cd907cb7647c" authors = ["Judith Schepers", "Benedikt V. Ehinger"] -version = "1.0.0-DEV" +version = "0.2.0" [deps] BSplineKit = "093aae92-e908-43d7-9660-e50ee39d5a0a" @@ -29,7 +29,7 @@ ClusterDepth = "0.2" MixedModels = "4" MixedModelsPermutations = "0.2" StatsModels = "0.7" -Unfold = "0.6,0.7" +Unfold = "0.7" julia = "1.9,1.10" [extras] From 3b4882cd4f68ea69c9007676f45b25657591bf55 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Sun, 23 Jun 2024 21:17:49 +0000 Subject: [PATCH 12/34] removed revise requirement --- test/setup.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/setup.jl b/test/setup.jl index 11ca32f..8a72add 100644 --- a/test/setup.jl +++ b/test/setup.jl @@ -4,5 +4,5 @@ using Unfold using DataFrames using StableRNGs -includet("../benchmark/types.jl") -includet("../benchmark/sim_and_fit.jl") +include("../benchmark/types.jl") +include("../benchmark/sim_and_fit.jl") From 1594e48a52a133db113835699326794b715f7b77 Mon Sep 17 00:00:00 2001 From: Benedikt Ehinger Date: Mon, 24 Jun 2024 09:04:09 +0200 Subject: [PATCH 13/34] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ccc7a70..403f4a0 100644 --- a/Project.toml +++ b/Project.toml @@ -29,7 +29,7 @@ ClusterDepth = "0.2" MixedModels = "4" MixedModelsPermutations = "0.2" StatsModels = "0.7" -Unfold = "0.7" +Unfold = "0.7.4" julia = "1.9,1.10" [extras] From 5d4bfbf421aeba3ca918297d7637e275039f639b Mon Sep 17 00:00:00 2001 From: Benedikt Ehinger Date: Mon, 24 Jun 2024 14:31:53 +0200 Subject: [PATCH 14/34] Update CI.yml --- .github/workflows/CI.yml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index ef00fe1..7744f1f 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -25,12 +25,12 @@ jobs: arch: - x64 steps: - - uses: actions/checkout@v3 - - uses: julia-actions/setup-julia@v1 + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - - uses: julia-actions/cache@v1 + - uses: julia-actions/cache@v2 - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 @@ -44,8 +44,8 @@ jobs: contents: write statuses: write steps: - - uses: actions/checkout@v3 - - uses: julia-actions/setup-julia@v1 + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 with: version: '1' - name: Configure doc environment From 096d0fa1a62af8fac8188475ff59203e9b1d8998 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Mon, 24 Jun 2024 13:21:36 +0000 Subject: [PATCH 15/34] fix docs mistake --- docs/make.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index b5c1186..58072b4 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -33,9 +33,10 @@ makedocs(; ), pages = [ "Home" => "index.md", - "Tutorials" => - ["Two-stage EEG analysis" => "generated/tutorial/two_stage_analysis.md"], - ["LMM EEG + ClusterDepth" => "generated/tutorial/lmm_clusterdepth.md"], + "Tutorials" => [ + "Two-stage EEG analysis" => "generated/tutorial/two_stage_analysis.md", + "LMM EEG + ClusterDepth" => "generated/tutorial/lmm_clusterdepth.md", + ], ], ) From d0fb7cc349c06bceac01c8071e3c8ee96b5f127d Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Mon, 24 Jun 2024 18:39:28 +0000 Subject: [PATCH 16/34] add statsmodels, update --- docs/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/Project.toml b/docs/Project.toml index 9c18e19..f880975 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -12,6 +12,7 @@ MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316" MixedModelsPermutations = "647c4018-d7ef-4d03-a0cc-8889a722319e" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" Unfold = "181c99d8-e21b-4ff3-b70b-c233eddec679" UnfoldMakie = "69a5ce3b-64fb-4f22-ae69-36dd4416af2a" UnfoldSim = "ed8ae6d2-84d3-44c6-ab46-0baf21700804" From 2589578f68ebafc0738431cf4a543d620e921140 Mon Sep 17 00:00:00 2001 From: Benedikt Ehinger Date: Fri, 10 Jan 2025 20:55:41 +0100 Subject: [PATCH 17/34] Update CI.yml --- .github/workflows/CI.yml | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 7744f1f..f8989c5 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -64,3 +64,9 @@ jobs: using UnfoldStats DocMeta.setdocmeta!(UnfoldStats, :DocTestSetup, :(using UnfoldStats); recursive=true) doctest(UnfoldStats)' + - name: Update UnfoldDocs + uses: peter-evans/repository-dispatch@v3 + with: + token: ${{ secrets.UNFOLDDOCSPAT }} + repository: unfoldtoolbox/UnfoldDocs + event-type: individual-docs-updated From d85d35d37c8e07b991d36e28d4b6fa15394fad77 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Fri, 10 Jan 2025 21:38:27 +0100 Subject: [PATCH 18/34] doc fixes --- docs/literate/tutorial/lmm_clusterdepth.jl | 1 + docs/literate/tutorial/two_stage_analysis.jl | 7 +- docs/src/EEGstats.md | 90 ++++++++++++++++++++ 3 files changed, 95 insertions(+), 3 deletions(-) create mode 100644 docs/src/EEGstats.md diff --git a/docs/literate/tutorial/lmm_clusterdepth.jl b/docs/literate/tutorial/lmm_clusterdepth.jl index 98d4e85..aec8c13 100644 --- a/docs/literate/tutorial/lmm_clusterdepth.jl +++ b/docs/literate/tutorial/lmm_clusterdepth.jl @@ -6,6 +6,7 @@ using MixedModelsPermutations, ClusterDepth # both necessary to activate correct using UnfoldStats using StatsModels using Random +using MixedModels srate = 25 design = MultiSubjectDesign(; diff --git a/docs/literate/tutorial/two_stage_analysis.jl b/docs/literate/tutorial/two_stage_analysis.jl index 7e28f52..8c17513 100644 --- a/docs/literate/tutorial/two_stage_analysis.jl +++ b/docs/literate/tutorial/two_stage_analysis.jl @@ -11,6 +11,7 @@ using HypothesisTests using CairoMakie using UnfoldMakie + # ## 1. Simulate data # This section can be skipped, if one already has (real) data that they want to analyse. @@ -124,7 +125,7 @@ effects_all_subjects = combine( ## Aggregate the marginal effects (per event type, time point and channel) over subjects ## using the mean as aggregation function aggregated_effects = @chain effects_all_subjects begin - groupby([:basisname, :channel, collect(keys(predictor_dict))..., :time]) + groupby([:eventname, :channel, collect(keys(predictor_dict))..., :time]) combine(:yhat .=> [x -> mean(skipmissing(x))] .=> Symbol("yhat_", mean)) end; first(aggregated_effects, 5) @@ -172,8 +173,8 @@ current_figure() ## Extract the spline coefficients for the `continuous` predictor variable ## from the Unfold models of all subjects -basisname = unique(effects_all_subjects.basisname) -coefs = extract_coefs(models.unfoldmodel, :continuous, basisname) +eventname = unique(effects_all_subjects.eventname) +coefs = extract_coefs(models.unfoldmodel, :continuous, eventname) ## Conduct a one-sample Hotelling's T² test separately for all channels and time points ## and compute a p-value. We compare the spline coefficients vector against 0. diff --git a/docs/src/EEGstats.md b/docs/src/EEGstats.md new file mode 100644 index 0000000..e90725b --- /dev/null +++ b/docs/src/EEGstats.md @@ -0,0 +1,90 @@ +````@example eeg-multichannel +using Random +using CairoMakie +using UnfoldSim +using Unfold +using UnfoldMakie +using Statistics +using ClusterDepth +```` + +# # Unfold single parameter testing + +This is an adaptation of the [ClusterDepth.jl tutorial](https://www.s-ccs.de/ClusterDepth.jl/dev/tutorials/eeg-multichannel/). + +Let's create data from 20 subjects +````@example eeg-multichannel +data,events = UnfoldSim.predef_eeg(20;return_epoched=true) +times = range(0,step=1/100,length=size(data,2)) +```` +Fit an UnfoldModel to each subject +````@example eeg +formula = @formula(0 ~ 1 + condition) +models = map((d, ev) -> (fit(UnfoldModel, formula, DataFrame(ev), d, ), ev.subject[1]), + eachslice(data; dims=3), + groupby(events, :subject)) +```` + +now we can inspect the data easily, and extract the face-effect + +````@example eeg +function add_subject!(df, s) + df[!, :subject] .= s + return df +end +allEffects = map((x) -> (effects(Dict(:condition => ["car", "face"]), x[1]), x[2]) |> (x) -> add_subject!(x[1], x[2]), models) |> e -> reduce(vcat, e) + +plot_erp(allEffects; mapping=(color=:condition, group=:subject)) +```` + +extract the face-coefficient from the linear model + +````@example eeg +allCoefs = map(m -> (coeftable(m[1]), m[2]) |> (x) -> add_subject!(x[1], x[2]), models) |> e -> reduce(vcat, e) +plot_erp(allCoefs; mapping=(group=:subject, col=:coefname)) +```` + +let's unstack the tidy-coef table into a matrix and put it to clusterdepth for clusterpermutation testing + +````@example eeg +faceCoefs = allCoefs |> x -> subset(x, :coefname => x -> x .== "condition: face") +erpMatrix = unstack(faceCoefs, :subject, :time, :estimate) |> x -> Matrix(x[:, 2:end])' |> collect +summary(erpMatrix) +```` + +## Clusterdepth + +````@example eeg +pvals = clusterdepth(erpMatrix; τ=quantile(TDist(n_subjects - 1), 0.95), nperm=5000); +nothing #hide +```` + +well - that was fast, less than a second for a cluster permutation test. not bad at all! + +## Plotting +Some plotting, and we add the identified cluster + +first calculate the ERP + +````@example eeg +faceERP = groupby(faceCoefs, [:time, :coefname]) |> + x -> combine(x, :estimate => mean => :estimate, + :estimate => std => :stderror); +nothing #hide +```` + +put the pvalues into a nicer format + +````@example eeg +pvalDF = ClusterDepth.cluster(pvals .<= 0.05) |> x -> DataFrame(:from => x[1] ./ 250, :to => (x[1] .+ x[2]) ./ 250, :coefname => "condition: face") +plot_erp(faceERP; stderror=true, pvalue=pvalDF) +```` + +Looks good to me! We identified the cluster :-) + +old unused code to use extra=(;pvalue=pvalDF) in the plotting function, but didnt work. + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + From 35b6e7d110dc6670f878e3ee04c4406a9306c268 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Thu, 23 Jan 2025 20:23:48 +0100 Subject: [PATCH 19/34] major update, might work with unfold0.8? --- Project.toml | 13 +++++----- docs/Project.toml | 2 +- ext/UnfoldStatsMixedModelsExt.jl | 13 ++++------ ext/UnfoldStatsMixedModelsPermutationsExt.jl | 25 +++++++++----------- src/UnfoldStats.jl | 3 ++- src/pvalue.jl | 3 +++ 6 files changed, 28 insertions(+), 31 deletions(-) create mode 100644 src/pvalue.jl diff --git a/Project.toml b/Project.toml index 403f4a0..91e846a 100644 --- a/Project.toml +++ b/Project.toml @@ -16,21 +16,22 @@ Unfold = "181c99d8-e21b-4ff3-b70b-c233eddec679" [weakdeps] ClusterDepth = "c8d8bbfa-f476-4995-adff-2987f04015d1" -MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316" MixedModelsPermutations = "647c4018-d7ef-4d03-a0cc-8889a722319e" +UnfoldMixedModels = "019ae9e0-8363-565c-86e5-97a5a2fe84f4" [extensions] -UnfoldStatsMixedModelsExt = ["MixedModels"] -UnfoldStatsMixedModelsPermutationsExt = ["ClusterDepth", "MixedModelsPermutations"] +UnfoldStatsMixedModelsExt = ["UnfoldMixedModels"] +UnfoldStatsMixedModelsPermutationsExt = ["UnfoldMixedModels", "ClusterDepth", "MixedModelsPermutations"] [compat] BSplineKit = "0.16,0.17" ClusterDepth = "0.2" -MixedModels = "4" MixedModelsPermutations = "0.2" +StatsAPI = "1.7.0" StatsModels = "0.7" -Unfold = "0.7.4" -julia = "1.9,1.10" +Unfold = "0.8" +UnfoldMixedModels = "0.1" +julia = "1.9,1.10,1.11" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/docs/Project.toml b/docs/Project.toml index f880975..c99c873 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -8,12 +8,12 @@ Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589" -MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316" MixedModelsPermutations = "647c4018-d7ef-4d03-a0cc-8889a722319e" StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" Unfold = "181c99d8-e21b-4ff3-b70b-c233eddec679" UnfoldMakie = "69a5ce3b-64fb-4f22-ae69-36dd4416af2a" +UnfoldMixedModels = "019ae9e0-8363-565c-86e5-97a5a2fe84f4" UnfoldSim = "ed8ae6d2-84d3-44c6-ab46-0baf21700804" UnfoldStats = "96fd419a-8306-4ce8-ba5b-cd907cb7647c" diff --git a/ext/UnfoldStatsMixedModelsExt.jl b/ext/UnfoldStatsMixedModelsExt.jl index 632373e..248a56d 100644 --- a/ext/UnfoldStatsMixedModelsExt.jl +++ b/ext/UnfoldStatsMixedModelsExt.jl @@ -1,23 +1,18 @@ module UnfoldStatsMixedModelsExt -using Unfold +using UnfoldMixedModels using UnfoldStats -using MixedModels -lmm_ext = Base.get_extension(Unfold, :UnfoldMixedModelsExt) - -if isnothing(lmm_ext) - error("Something went wrong with getting the Unfold UnfoldMixedModelsExt extension") -end # Currently, `extract_coefs` is not implemented for mixed-effects models UnfoldStats.extract_coefs( model::Union{ - lmm_ext.UnfoldLinearMixedModel, - lmm_ext.UnfoldLinearMixedModelContinuousTime, + UnfoldLinearMixedModel, + UnfoldLinearMixedModelContinuousTime, }, predictor, basisname, ) = throw( "The `extract_coefs` function is currently not implemented for mixed-effects models.", ) + end \ No newline at end of file diff --git a/ext/UnfoldStatsMixedModelsPermutationsExt.jl b/ext/UnfoldStatsMixedModelsPermutationsExt.jl index 2132b67..be40c16 100644 --- a/ext/UnfoldStatsMixedModelsPermutationsExt.jl +++ b/ext/UnfoldStatsMixedModelsPermutationsExt.jl @@ -1,19 +1,16 @@ module UnfoldStatsMixedModelsPermutationsExt -using Unfold -import Unfold: pvalues + +using UnfoldMixedModels using UnfoldStats #using MixedModels using MixedModelsPermutations using ClusterDepth using Logging using Random -const MixedModels = MixedModelsPermutations.MixedModels -LMMext = Base.get_extension(Unfold, :UnfoldMixedModelsExt) - -function Unfold.pvalues( +function UnfoldStats.pvalue( rng, - model::LMMext.UnfoldLinearMixedModel, + model::UnfoldLinearMixedModel, data, coefficient; type = "clusterdepth", @@ -99,9 +96,13 @@ function lmm_permutations( permdata = Matrix{Float64}(undef, length(time_selection), n_permutations) - Xs = LMMext.prepare_modelmatrix(model) + Xs = UnfoldMixedModels.prepare_modelmatrix(model) - mm_outer = LMMext.LinearMixedModel_wrapper(Unfold.formulas(model), data[1, 1, :], Xs) + mm_outer = UnfoldMixedModels.LinearMixedModel_wrapper( + Unfold.formulas(model), + data[1, 1, :], + Xs, + ) mm_outer.optsum.maxtime = 0.1 # chIx = 1 # for now we only support 1 channel anyway @@ -160,11 +161,7 @@ function get_lmm_statistic( ] end -function get_lmm_statistic( - model::LMMext.UnfoldLinearMixedModel, - coefficient::Int, - lmm_statistic, -) +function get_lmm_statistic(model::UnfoldLinearMixedModel, coefficient::Int, lmm_statistic) return get_lmm_statistic(model, modelfit(model), coefficient, lmm_statistic) # r = coeftable(m) # r = subset(r, :group => (x -> isnothing.(x)), :coefname => (x -> x .!== "(Intercept)")) diff --git a/src/UnfoldStats.jl b/src/UnfoldStats.jl index 540dcf5..d8f0140 100644 --- a/src/UnfoldStats.jl +++ b/src/UnfoldStats.jl @@ -4,7 +4,8 @@ using BSplineKit: Reexport using Unfold using BSplineKit using StatsModels - +using StatsAPI +include("pvalue.jl") include("extract_coefs.jl") # export functions to extract model coefficients diff --git a/src/pvalue.jl b/src/pvalue.jl new file mode 100644 index 0000000..ace64cc --- /dev/null +++ b/src/pvalue.jl @@ -0,0 +1,3 @@ +# needed for pvalues in the extension +@deprecate pvalues(args...) pvalue(args...) +pvalue(args...; kwargs...) = error("not implemented, see `methods(pvalues)` for a list") \ No newline at end of file From 6d1f6727bb119829c7b09589848a1422349b92a1 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Thu, 23 Jan 2025 20:24:29 +0100 Subject: [PATCH 20/34] compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 91e846a..3bd06b4 100644 --- a/Project.toml +++ b/Project.toml @@ -27,7 +27,7 @@ UnfoldStatsMixedModelsPermutationsExt = ["UnfoldMixedModels", "ClusterDepth", "M BSplineKit = "0.16,0.17" ClusterDepth = "0.2" MixedModelsPermutations = "0.2" -StatsAPI = "1.7.0" +StatsAPI = "1" StatsModels = "0.7" Unfold = "0.8" UnfoldMixedModels = "0.1" From b58a866c2aa9818a7cf98f3c9e26ea45af108dab Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Fri, 24 Jan 2025 00:04:28 +0100 Subject: [PATCH 21/34] fix a lot of bugs for newer unfold versions --- benchmark/sim_and_fit.jl | 6 ++--- docs/literate/tutorial/lmm_clusterdepth.jl | 6 ++--- docs/literate/tutorial/two_stage_analysis.jl | 26 ++++++++++---------- src/UnfoldStats.jl | 2 ++ src/pvalue.jl | 5 ++-- 5 files changed, 24 insertions(+), 21 deletions(-) diff --git a/benchmark/sim_and_fit.jl b/benchmark/sim_and_fit.jl index 9a7a9db..aa38224 100644 --- a/benchmark/sim_and_fit.jl +++ b/benchmark/sim_and_fit.jl @@ -1,5 +1,5 @@ using UnfoldSim -using Unfold +using Unfold, UnfoldMixedModels using StableRNGs using DataFrames @@ -36,8 +36,8 @@ function sim_and_fit( ) # At the moment, the function is not implemented for mixed models - ext = Base.get_extension(Unfold, :UnfoldMixedModelsExt) - if model_type in [ext.UnfoldLinearMixedModel, ext.UnfoldLinearMixedModelContinuousTime] + #ext = Base.get_extension(Unfold, :UnfoldMixedModelsExt) + if model_type in [UnfoldLinearMixedModel, UnfoldLinearMixedModelContinuousTime] throw("Not implemented.") end diff --git a/docs/literate/tutorial/lmm_clusterdepth.jl b/docs/literate/tutorial/lmm_clusterdepth.jl index aec8c13..6b1c577 100644 --- a/docs/literate/tutorial/lmm_clusterdepth.jl +++ b/docs/literate/tutorial/lmm_clusterdepth.jl @@ -1,12 +1,12 @@ # get some data using UnfoldSim -using Unfold +using UnfoldMixedModels + using MixedModelsPermutations, ClusterDepth # both necessary to activate correct extension! using UnfoldStats using StatsModels using Random -using MixedModels srate = 25 design = MultiSubjectDesign(; @@ -71,7 +71,7 @@ m = fit( # # Cluster Permute :) coefficient = 2 -pvalues( +pvalue( MersenneTwister(1), m, data_e, diff --git a/docs/literate/tutorial/two_stage_analysis.jl b/docs/literate/tutorial/two_stage_analysis.jl index 8c17513..ea9e97a 100644 --- a/docs/literate/tutorial/two_stage_analysis.jl +++ b/docs/literate/tutorial/two_stage_analysis.jl @@ -56,7 +56,7 @@ noise = PinkNoise(; noiselevel = 5) # # ``` -data, events = simulate( +data, events = UnfoldSim.simulate( StableRNG(1), design, signal_multichannel, @@ -84,7 +84,7 @@ basisfunction = firbasis((-0.1, 0.7), 100) formula = @formula 0 ~ 1 + spl(continuous, 4) ## Combine basisfunction and formula in an event dict -event_dict = Dict(Any => (formula, basisfunction)) +event_vec = [Any => (formula, basisfunction)] subject_list = unique(events.subject) model_list = UnfoldLinearModelContinuousTime[] @@ -95,7 +95,7 @@ data_slices = eachslice(data, dims = ndims(data)) for s = 1:size(data, ndims(data)) m = fit( UnfoldModel, - event_dict, + event_vec, subset(events, :subject => ByRow(==(subject_list[s]))), data_slices[s], ) @@ -144,15 +144,15 @@ pos2d = [Point2f(p[1] + 0.5, p[2] + 0.5) for p in pos2d]; # The rows (from top to bottom) represent the marginal effects for different levels of the predictor `continuous`. ## Set the size of the time bins for the topoplot series -bin_size = 0.1 - +bin_width = 0.1 +aggregated_effects.estimate = aggregated_effects.yhat_mean # bug in UnfoldMakie 0.5.12 f_effects = Figure(size = (1200, 600)) tp_effects = plot_topoplotseries!( f_effects, - aggregated_effects, - bin_size, + aggregated_effects; + bin_width, positions = pos2d, - mapping = (; y = :yhat_mean, row = :continuous), + mapping = (; y = :estimate, row = :continuous), visual = (; enlarge = 0.6, label_scatter = false, colorrange = (-3, 3)), ) @@ -195,15 +195,15 @@ p_values_df = DataFrame( first(p_values_df, 5) # As a last step, we visualize the p-values in a topoplot series. -bin_size = 0.1 - +bin_width = 0.1 +p_values_df.estimate = p_values_df.p_values f_pvalues = Figure(size = (1200, 200)) tp_pvalues = plot_topoplotseries!( f_pvalues, - p_values_df, - bin_size, + p_values_df; + bin_width, positions = pos2d, - mapping = (; y = :p_values), + mapping = (; y = :estimate), visual = (; enlarge = 0.6, label_scatter = false, colorrange = (0, 0.1)), colorbar = (; label = "p-value"), ) diff --git a/src/UnfoldStats.jl b/src/UnfoldStats.jl index d8f0140..0b7f8de 100644 --- a/src/UnfoldStats.jl +++ b/src/UnfoldStats.jl @@ -5,10 +5,12 @@ using Unfold using BSplineKit using StatsModels using StatsAPI +import StatsAPI: pvalue include("pvalue.jl") include("extract_coefs.jl") # export functions to extract model coefficients export extract_coefs +export pvalue end diff --git a/src/pvalue.jl b/src/pvalue.jl index ace64cc..885a196 100644 --- a/src/pvalue.jl +++ b/src/pvalue.jl @@ -1,3 +1,4 @@ # needed for pvalues in the extension -@deprecate pvalues(args...) pvalue(args...) -pvalue(args...; kwargs...) = error("not implemented, see `methods(pvalues)` for a list") \ No newline at end of file +#@deprecate pvalues(args...) pvalue(args...) +StatsAPI.pvalue(args...; kwargs...) = + error("not implemented, see `methods(pvalues)` for a list") \ No newline at end of file From d2670dafc1c0c393765adec1c54c2c8c6c3473a8 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Fri, 24 Jan 2025 09:49:08 +0100 Subject: [PATCH 22/34] fix fix fix everyhtinggit add test --- ...{two_stage_analysis.jl => test_splines.jl} | 0 docs/make.jl | 5 +++-- .../test_single_coefficient.md} | 21 +++++++------------ test/Project.toml | 1 + 4 files changed, 11 insertions(+), 16 deletions(-) rename docs/literate/tutorial/{two_stage_analysis.jl => test_splines.jl} (100%) rename docs/src/{EEGstats.md => tutorial/test_single_coefficient.md} (87%) diff --git a/docs/literate/tutorial/two_stage_analysis.jl b/docs/literate/tutorial/test_splines.jl similarity index 100% rename from docs/literate/tutorial/two_stage_analysis.jl rename to docs/literate/tutorial/test_splines.jl diff --git a/docs/make.jl b/docs/make.jl index 58072b4..eae7d1e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -34,8 +34,9 @@ makedocs(; pages = [ "Home" => "index.md", "Tutorials" => [ - "Two-stage EEG analysis" => "generated/tutorial/two_stage_analysis.md", - "LMM EEG + ClusterDepth" => "generated/tutorial/lmm_clusterdepth.md", + "2-stage single parameter test" => "tutorial/test_single_coefficient.md", + "Test multi-channel spline-effects" => "generated/tutorial/two_stage_analysis.md", + "**MixedModels** EEG + Clusterpermutation" => "generated/tutorial/lmm_clusterdepth.md", ], ], ) diff --git a/docs/src/EEGstats.md b/docs/src/tutorial/test_single_coefficient.md similarity index 87% rename from docs/src/EEGstats.md rename to docs/src/tutorial/test_single_coefficient.md index e90725b..857fb25 100644 --- a/docs/src/EEGstats.md +++ b/docs/src/tutorial/test_single_coefficient.md @@ -27,7 +27,7 @@ models = map((d, ev) -> (fit(UnfoldModel, formula, DataFrame(ev), d, ), ev.subje now we can inspect the data easily, and extract the face-effect -````@example eeg +````@example eeg-multichannel function add_subject!(df, s) df[!, :subject] .= s return df @@ -39,14 +39,14 @@ plot_erp(allEffects; mapping=(color=:condition, group=:subject)) extract the face-coefficient from the linear model -````@example eeg +````@example eeg-multichannel allCoefs = map(m -> (coeftable(m[1]), m[2]) |> (x) -> add_subject!(x[1], x[2]), models) |> e -> reduce(vcat, e) plot_erp(allCoefs; mapping=(group=:subject, col=:coefname)) ```` let's unstack the tidy-coef table into a matrix and put it to clusterdepth for clusterpermutation testing -````@example eeg +````@example eeg-multichannel faceCoefs = allCoefs |> x -> subset(x, :coefname => x -> x .== "condition: face") erpMatrix = unstack(faceCoefs, :subject, :time, :estimate) |> x -> Matrix(x[:, 2:end])' |> collect summary(erpMatrix) @@ -54,7 +54,7 @@ summary(erpMatrix) ## Clusterdepth -````@example eeg +````@example eeg-multichannel pvals = clusterdepth(erpMatrix; τ=quantile(TDist(n_subjects - 1), 0.95), nperm=5000); nothing #hide ```` @@ -66,7 +66,7 @@ Some plotting, and we add the identified cluster first calculate the ERP -````@example eeg +````@example eeg-multichannel faceERP = groupby(faceCoefs, [:time, :coefname]) |> x -> combine(x, :estimate => mean => :estimate, :estimate => std => :stderror); @@ -75,16 +75,9 @@ nothing #hide put the pvalues into a nicer format -````@example eeg +````@example eeg-multichannel pvalDF = ClusterDepth.cluster(pvals .<= 0.05) |> x -> DataFrame(:from => x[1] ./ 250, :to => (x[1] .+ x[2]) ./ 250, :coefname => "condition: face") plot_erp(faceERP; stderror=true, pvalue=pvalDF) ```` -Looks good to me! We identified the cluster :-) - -old unused code to use extra=(;pvalue=pvalDF) in the plotting function, but didnt work. - ---- - -*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* - +Looks good to me! We identified the cluster :-) \ No newline at end of file diff --git a/test/Project.toml b/test/Project.toml index 78d6ee2..12a460f 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -4,4 +4,5 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" StatsModels = "3eaba693-59b7-5ba5-a881-562e759f1c8d" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Unfold = "181c99d8-e21b-4ff3-b70b-c233eddec679" +UnfoldMixedModels = "019ae9e0-8363-565c-86e5-97a5a2fe84f4" UnfoldSim = "ed8ae6d2-84d3-44c6-ab46-0baf21700804" From 1bec33bcc63c67eb92af67f644e557c153b68965 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Fri, 24 Jan 2025 10:12:57 +0100 Subject: [PATCH 23/34] forgot to rename --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index eae7d1e..ddb3b19 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -35,7 +35,7 @@ makedocs(; "Home" => "index.md", "Tutorials" => [ "2-stage single parameter test" => "tutorial/test_single_coefficient.md", - "Test multi-channel spline-effects" => "generated/tutorial/two_stage_analysis.md", + "Test multi-channel spline-effects" => "generated/tutorial/test_splines.md", "**MixedModels** EEG + Clusterpermutation" => "generated/tutorial/lmm_clusterdepth.md", ], ], From dd44fada6489cefbe775273eaf7c4e3d1c2087d5 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Fri, 24 Jan 2025 21:24:32 +0100 Subject: [PATCH 24/34] some conflicts with MixedModels.simulate and UnfoldSim.simulate --- benchmark/sim_and_fit.jl | 6 ++++-- docs/src/tutorial/test_single_coefficient.md | 18 +++++++++--------- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/benchmark/sim_and_fit.jl b/benchmark/sim_and_fit.jl index aa38224..1bf0f0a 100644 --- a/benchmark/sim_and_fit.jl +++ b/benchmark/sim_and_fit.jl @@ -82,7 +82,8 @@ simulate_data(sim_type::T, simulation, return_epoched, seed) where {T} = function simulate_data(::SingleEventType, simulation, return_epoched, seed) # Simulate data - data, events = simulate(StableRNG(seed), simulation; return_epoched = return_epoched) + data, events = + UnfoldSim.simulate(StableRNG(seed), simulation; return_epoched = return_epoched) return data, events end @@ -90,7 +91,8 @@ end function simulate_data(::MultipleEventTypes, simulation, return_epoched, seed) # Simulate data - data, events = simulate(StableRNG(seed), simulation; return_epoched = return_epoched) + data, events = + UnfoldSim.simulate(StableRNG(seed), simulation; return_epoched = return_epoched) # Add an event column to the events df and assign each event to an event type events[!, :event] = repeat(["stim", "fix"], size(events, 1) ÷ 2) diff --git a/docs/src/tutorial/test_single_coefficient.md b/docs/src/tutorial/test_single_coefficient.md index 857fb25..bbc3955 100644 --- a/docs/src/tutorial/test_single_coefficient.md +++ b/docs/src/tutorial/test_single_coefficient.md @@ -1,4 +1,4 @@ -````@example eeg-multichannel +````@example eeg using Random using CairoMakie using UnfoldSim @@ -10,10 +10,10 @@ using ClusterDepth # # Unfold single parameter testing -This is an adaptation of the [ClusterDepth.jl tutorial](https://www.s-ccs.de/ClusterDepth.jl/dev/tutorials/eeg-multichannel/). +This is an adaptation of the [ClusterDepth.jl tutorial](https://www.s-ccs.de/ClusterDepth.jl/dev/tutorials/eeg/). Let's create data from 20 subjects -````@example eeg-multichannel +````@example eeg data,events = UnfoldSim.predef_eeg(20;return_epoched=true) times = range(0,step=1/100,length=size(data,2)) ```` @@ -27,7 +27,7 @@ models = map((d, ev) -> (fit(UnfoldModel, formula, DataFrame(ev), d, ), ev.subje now we can inspect the data easily, and extract the face-effect -````@example eeg-multichannel +````@example eeg function add_subject!(df, s) df[!, :subject] .= s return df @@ -39,14 +39,14 @@ plot_erp(allEffects; mapping=(color=:condition, group=:subject)) extract the face-coefficient from the linear model -````@example eeg-multichannel +````@example eeg allCoefs = map(m -> (coeftable(m[1]), m[2]) |> (x) -> add_subject!(x[1], x[2]), models) |> e -> reduce(vcat, e) plot_erp(allCoefs; mapping=(group=:subject, col=:coefname)) ```` let's unstack the tidy-coef table into a matrix and put it to clusterdepth for clusterpermutation testing -````@example eeg-multichannel +````@example eeg faceCoefs = allCoefs |> x -> subset(x, :coefname => x -> x .== "condition: face") erpMatrix = unstack(faceCoefs, :subject, :time, :estimate) |> x -> Matrix(x[:, 2:end])' |> collect summary(erpMatrix) @@ -54,7 +54,7 @@ summary(erpMatrix) ## Clusterdepth -````@example eeg-multichannel +````@example eeg pvals = clusterdepth(erpMatrix; τ=quantile(TDist(n_subjects - 1), 0.95), nperm=5000); nothing #hide ```` @@ -66,7 +66,7 @@ Some plotting, and we add the identified cluster first calculate the ERP -````@example eeg-multichannel +````@example eeg faceERP = groupby(faceCoefs, [:time, :coefname]) |> x -> combine(x, :estimate => mean => :estimate, :estimate => std => :stderror); @@ -75,7 +75,7 @@ nothing #hide put the pvalues into a nicer format -````@example eeg-multichannel +````@example eeg pvalDF = ClusterDepth.cluster(pvals .<= 0.05) |> x -> DataFrame(:from => x[1] ./ 250, :to => (x[1] .+ x[2]) ./ 250, :coefname => "condition: face") plot_erp(faceERP; stderror=true, pvalue=pvalDF) ```` From 4e3331d5128ca7cf79aae9202a855870e63197e4 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Mon, 27 Jan 2025 11:52:57 +0100 Subject: [PATCH 25/34] finally found the bug --- docs/src/tutorial/test_single_coefficient.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/docs/src/tutorial/test_single_coefficient.md b/docs/src/tutorial/test_single_coefficient.md index bbc3955..74df36a 100644 --- a/docs/src/tutorial/test_single_coefficient.md +++ b/docs/src/tutorial/test_single_coefficient.md @@ -6,6 +6,7 @@ using Unfold using UnfoldMakie using Statistics using ClusterDepth +using DataFrames ```` # # Unfold single parameter testing @@ -15,12 +16,12 @@ This is an adaptation of the [ClusterDepth.jl tutorial](https://www.s-ccs.de/Clu Let's create data from 20 subjects ````@example eeg data,events = UnfoldSim.predef_eeg(20;return_epoched=true) -times = range(0,step=1/100,length=size(data,2)) +times = range(0,step=1/100,length=size(data,1)) ```` Fit an UnfoldModel to each subject ````@example eeg formula = @formula(0 ~ 1 + condition) -models = map((d, ev) -> (fit(UnfoldModel, formula, DataFrame(ev), d, ), ev.subject[1]), +models = map((d, ev) -> (fit(UnfoldModel, formula, DataFrame(ev), collect(d), times), ev.subject[1]), eachslice(data; dims=3), groupby(events, :subject)) ```` From baf213b93dc23443287a4c3ebaffd94ea966a4fe Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Mon, 27 Jan 2025 17:46:49 +0100 Subject: [PATCH 26/34] maybe maybe maybe --- docs/Project.toml | 1 + docs/src/tutorial/test_single_coefficient.md | 1 + 2 files changed, 2 insertions(+) diff --git a/docs/Project.toml b/docs/Project.toml index c99c873..b8f6a83 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -3,6 +3,7 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Chain = "8be319e6-bccf-4806-a6f7-6fae938471bc" ClusterDepth = "c8d8bbfa-f476-4995-adff-2987f04015d1" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5" diff --git a/docs/src/tutorial/test_single_coefficient.md b/docs/src/tutorial/test_single_coefficient.md index 74df36a..1dbe993 100644 --- a/docs/src/tutorial/test_single_coefficient.md +++ b/docs/src/tutorial/test_single_coefficient.md @@ -7,6 +7,7 @@ using UnfoldMakie using Statistics using ClusterDepth using DataFrames +using Distributions ```` # # Unfold single parameter testing From 0f23a859339f03cbe36333d4c9bd7f506da29569 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Mon, 27 Jan 2025 17:47:16 +0100 Subject: [PATCH 27/34] mayyybe --- docs/src/tutorial/test_single_coefficient.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/docs/src/tutorial/test_single_coefficient.md b/docs/src/tutorial/test_single_coefficient.md index 1dbe993..d71a9f6 100644 --- a/docs/src/tutorial/test_single_coefficient.md +++ b/docs/src/tutorial/test_single_coefficient.md @@ -16,7 +16,8 @@ This is an adaptation of the [ClusterDepth.jl tutorial](https://www.s-ccs.de/Clu Let's create data from 20 subjects ````@example eeg -data,events = UnfoldSim.predef_eeg(20;return_epoched=true) +n_subjects = 20 +data,events = UnfoldSim.predef_eeg(n_subjects;return_epoched=true) times = range(0,step=1/100,length=size(data,1)) ```` Fit an UnfoldModel to each subject @@ -57,6 +58,7 @@ summary(erpMatrix) ## Clusterdepth ````@example eeg + pvals = clusterdepth(erpMatrix; τ=quantile(TDist(n_subjects - 1), 0.95), nperm=5000); nothing #hide ```` @@ -79,7 +81,7 @@ put the pvalues into a nicer format ````@example eeg pvalDF = ClusterDepth.cluster(pvals .<= 0.05) |> x -> DataFrame(:from => x[1] ./ 250, :to => (x[1] .+ x[2]) ./ 250, :coefname => "condition: face") -plot_erp(faceERP; stderror=true, pvalue=pvalDF) +plot_erp(faceERP; stderror=true, significance=pvalDF) ```` Looks good to me! We identified the cluster :-) \ No newline at end of file From 56746f5fd014e4bf19e1a448a1a161c4d1aa1001 Mon Sep 17 00:00:00 2001 From: Benedikt Ehinger Date: Thu, 13 Feb 2025 15:19:13 +0100 Subject: [PATCH 28/34] Update sim_and_fit.jl --- benchmark/sim_and_fit.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/benchmark/sim_and_fit.jl b/benchmark/sim_and_fit.jl index 1bf0f0a..9d89645 100644 --- a/benchmark/sim_and_fit.jl +++ b/benchmark/sim_and_fit.jl @@ -52,7 +52,7 @@ function sim_and_fit( data, events = simulate_data(sim_type, simulation, return_epoched, seed) # Create event dict containing basis function(s)/times and formula(s) for all events - event_dict = create_event_dict(sim_type, model_type, simulation) + event_dict = create_event_vector(sim_type, model_type, simulation) # Fit an Unfold model for each subject subject_list = unique(events.subject) @@ -100,14 +100,14 @@ function simulate_data(::MultipleEventTypes, simulation, return_epoched, seed) return data, events end -create_event_dict(sim_type::T, model_type, simulation) where {T} = create_event_dict( +create_event_vector(sim_type::T, model_type, simulation) where {T} = create_event_vector( EventStyle(T), PredictorStyle(T), model_type::Type{<:UnfoldModel}, simulation, ) -function create_event_dict( +function create_event_vector( ::MultipleEventTypes, ::ManyPredictors, model_type::Type{<:UnfoldModel}, @@ -134,7 +134,7 @@ function create_event_dict( return event_dict end -function create_event_dict( +function create_event_vector( ::MultipleEventTypes, ::OnlySplines, model_type::Type{<:UnfoldModel}, @@ -161,7 +161,7 @@ function create_event_dict( return event_vec end -function create_event_dict( +function create_event_vector( ::SingleEventType, ::ManyPredictors, model_type::Type{<:UnfoldModel}, @@ -184,7 +184,7 @@ function create_event_dict( return event_vec end -function create_event_dict( +function create_event_vector( ::SingleEventType, ::OnlySplines, model_type::Type{<:UnfoldModel}, From e7a0a13f8e015c259bed8b4b05730a1fb144a44b Mon Sep 17 00:00:00 2001 From: Benedikt Ehinger Date: Thu, 13 Feb 2025 15:42:47 +0100 Subject: [PATCH 29/34] Apply suggestions from code review @jschepers did it Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- benchmark/sim_and_fit.jl | 1 - docs/literate/tutorial/lmm_clusterdepth.jl | 3 ++- docs/make.jl | 4 ++-- ext/UnfoldStatsMixedModelsExt.jl | 7 ++----- src/UnfoldStats.jl | 1 - src/pvalue.jl | 3 +-- test/extract_coefs.jl | 3 +-- 7 files changed, 8 insertions(+), 14 deletions(-) diff --git a/benchmark/sim_and_fit.jl b/benchmark/sim_and_fit.jl index 9d89645..6ba43e1 100644 --- a/benchmark/sim_and_fit.jl +++ b/benchmark/sim_and_fit.jl @@ -36,7 +36,6 @@ function sim_and_fit( ) # At the moment, the function is not implemented for mixed models - #ext = Base.get_extension(Unfold, :UnfoldMixedModelsExt) if model_type in [UnfoldLinearMixedModel, UnfoldLinearMixedModelContinuousTime] throw("Not implemented.") end diff --git a/docs/literate/tutorial/lmm_clusterdepth.jl b/docs/literate/tutorial/lmm_clusterdepth.jl index 6b1c577..df47df1 100644 --- a/docs/literate/tutorial/lmm_clusterdepth.jl +++ b/docs/literate/tutorial/lmm_clusterdepth.jl @@ -1,4 +1,5 @@ - +# !!! important +# This functionality is not tested. While it returns reasonable results, we haven't written any unit-tests, nor tested the type-1 error probability yet # get some data using UnfoldSim using UnfoldMixedModels diff --git a/docs/make.jl b/docs/make.jl index ddb3b19..c9cb49f 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -34,8 +34,8 @@ makedocs(; pages = [ "Home" => "index.md", "Tutorials" => [ - "2-stage single parameter test" => "tutorial/test_single_coefficient.md", - "Test multi-channel spline-effects" => "generated/tutorial/test_splines.md", + "two stage single parameter test (ttest)" => "tutorial/test_single_coefficient.md", + "two stage multi parameter test (hotelling-T²)) => "generated/tutorial/test_splines.md", "**MixedModels** EEG + Clusterpermutation" => "generated/tutorial/lmm_clusterdepth.md", ], ], diff --git a/ext/UnfoldStatsMixedModelsExt.jl b/ext/UnfoldStatsMixedModelsExt.jl index 248a56d..f6630f4 100644 --- a/ext/UnfoldStatsMixedModelsExt.jl +++ b/ext/UnfoldStatsMixedModelsExt.jl @@ -5,14 +5,11 @@ using UnfoldStats # Currently, `extract_coefs` is not implemented for mixed-effects models UnfoldStats.extract_coefs( - model::Union{ - UnfoldLinearMixedModel, - UnfoldLinearMixedModelContinuousTime, - }, + model::Union{UnfoldLinearMixedModel,UnfoldLinearMixedModelContinuousTime}, predictor, basisname, ) = throw( "The `extract_coefs` function is currently not implemented for mixed-effects models.", ) - + end \ No newline at end of file diff --git a/src/UnfoldStats.jl b/src/UnfoldStats.jl index 0b7f8de..4030a35 100644 --- a/src/UnfoldStats.jl +++ b/src/UnfoldStats.jl @@ -1,6 +1,5 @@ module UnfoldStats -using BSplineKit: Reexport using Unfold using BSplineKit using StatsModels diff --git a/src/pvalue.jl b/src/pvalue.jl index 885a196..eac6bff 100644 --- a/src/pvalue.jl +++ b/src/pvalue.jl @@ -1,4 +1,3 @@ # needed for pvalues in the extension -#@deprecate pvalues(args...) pvalue(args...) StatsAPI.pvalue(args...; kwargs...) = - error("not implemented, see `methods(pvalues)` for a list") \ No newline at end of file + error("not implemented, see `methods(pvalue)` for a list") \ No newline at end of file diff --git a/test/extract_coefs.jl b/test/extract_coefs.jl index 7c855c6..63fac54 100644 --- a/test/extract_coefs.jl +++ b/test/extract_coefs.jl @@ -8,7 +8,6 @@ σs = Dict(:subject => [1, 0.1, 0.1, 0.1]) simulation = define_simulation(sim_type, β, σs, n_subjects = 5) - model_type = UnfoldLinearModelContinuousTime for model_type in [UnfoldLinearModel, UnfoldLinearModelContinuousTime] models = sim_and_fit(sim_type, simulation, model_type, seed = 1) @@ -44,7 +43,7 @@ # Test exceptions @test_throws ArgumentError extract_coefs(model_1, "continuous", "fix") - @test_throws ArgumentError extract_coefs(model_1, :continuous, "bla fix") + @test_throws ArgumentError extract_coefs(model_1, :continuous, "this-basis-doesnt-exist") @test_throws ArgumentError extract_coefs(model_1, :pet, "stim") # Test `extract_coefs` method for single Unfold model From 8ae7a9389905a5972134a6b086f37ee2d163df10 Mon Sep 17 00:00:00 2001 From: Judith Schepers Date: Thu, 13 Feb 2025 16:31:04 +0100 Subject: [PATCH 30/34] Add docs preview for PRs --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index c9cb49f..01b0b1a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -41,4 +41,4 @@ makedocs(; ], ) -deploydocs(; repo = "github.com/unfoldtoolbox/UnfoldStats.jl", devbranch = "main") +deploydocs(; repo = "github.com/unfoldtoolbox/UnfoldStats.jl", devbranch = "main", push_preview = true) From c12dd9193b3091e2db1a98d1841d85f233a3f6e4 Mon Sep 17 00:00:00 2001 From: Judith Schepers Date: Thu, 13 Feb 2025 17:02:25 +0100 Subject: [PATCH 31/34] Fix parantheses and capitalization --- docs/make.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 01b0b1a..bbf6008 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -34,9 +34,9 @@ makedocs(; pages = [ "Home" => "index.md", "Tutorials" => [ - "two stage single parameter test (ttest)" => "tutorial/test_single_coefficient.md", - "two stage multi parameter test (hotelling-T²)) => "generated/tutorial/test_splines.md", - "**MixedModels** EEG + Clusterpermutation" => "generated/tutorial/lmm_clusterdepth.md", + "Two-stage single parameter test (t-test)" => "tutorial/test_single_coefficient.md", + "Two-stage multi parameter test (Hotelling's T-squared test)" => "generated/tutorial/test_splines.md", + "**MixedModels** EEG + Cluster permutation test" => "generated/tutorial/lmm_clusterdepth.md", ], ], ) From 6b9ba7e2b3f627b36d7a23909059ae34bbbfb22d Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Tue, 18 Feb 2025 23:12:51 +0100 Subject: [PATCH 32/34] improved tutorials --- docs/literate/tutorial/lmm_clusterdepth.jl | 28 ++++++-- docs/src/tutorial/test_single_coefficient.md | 73 +++++++++++++------- 2 files changed, 71 insertions(+), 30 deletions(-) diff --git a/docs/literate/tutorial/lmm_clusterdepth.jl b/docs/literate/tutorial/lmm_clusterdepth.jl index df47df1..c60695d 100644 --- a/docs/literate/tutorial/lmm_clusterdepth.jl +++ b/docs/literate/tutorial/lmm_clusterdepth.jl @@ -1,6 +1,14 @@ # !!! important # This functionality is not tested. While it returns reasonable results, we haven't written any unit-tests, nor tested the type-1 error probability yet -# get some data + + +# ## 1. Simulate data +# This section can be skipped, if one already has (real) data that they want to analyse. + +# ```@raw html +#
+# Click to expand the simulation details +# ``` using UnfoldSim using UnfoldMixedModels @@ -53,10 +61,12 @@ data_e, events = UnfoldSim.simulate( times = range(-0.1, 0.5, length = size(data_e, 1)) data_e = reshape(data_e, 1, size(data_e, 1), :) #events.latency .+= repeat(range(0,length=size(data,2),step=size(data,1)),inner=size(events[events.subject.=="S01",:],1)) +# ```@raw html +#
+# ``` - - -# # Fit LMM +# ## 2. Fit mass-univariate LMMs +# We have some typical experimental data with subject and item effects. Item refer to stimuli here, based on our `stimtype` condition these are either different `cars` or `faces`. m = fit( UnfoldModel, [ @@ -70,8 +80,16 @@ m = fit( ); -# # Cluster Permute :) +# ## 3. Cluster permutation test +# If we would run a statistical test on each time-point separately, we would greatly inflate the type-1 error, reaching significance on any each sample much higher than the assumed α=0.05. +# One solution are cluster permutation test, where we instead test for clustersizes of connected significant clusters. In "classical" two-stage testing, such a permutation test is straight forward. But for LMMs we have to think of something more clever, as it is not directly clear how to permute if both subject and item effects exist (you gonna break the relation between the two). We did that in `MixedModelsPermutations` and can apply this strategy to EEG data as well.` + +# select the fixed-effects coefficient to test (`stimtype`) coefficient = 2 + +# call the permutatio test +# !!! note +# This interface is very likely to change in the future pvalue( MersenneTwister(1), m, diff --git a/docs/src/tutorial/test_single_coefficient.md b/docs/src/tutorial/test_single_coefficient.md index d71a9f6..14e197a 100644 --- a/docs/src/tutorial/test_single_coefficient.md +++ b/docs/src/tutorial/test_single_coefficient.md @@ -1,4 +1,18 @@ -````@example eeg + + +# Unfold single parameter cluster permutation testing +This is an adaptation of the [ClusterDepth.jl tutorial](https://www.s-ccs.de/ClusterDepth.jl/dev/tutorials/eeg/). + +We simulate a 1x2 design and perform cluster permutation testing via the `ClusterDepth.jl` package + +## Simulate data +Let's create data from 20 subjects + + ```@raw html +
+ Click to expand the simulation details + ``` +```@example eeg using Random using CairoMakie using UnfoldSim @@ -8,29 +22,34 @@ using Statistics using ClusterDepth using DataFrames using Distributions -```` +``` -# # Unfold single parameter testing -This is an adaptation of the [ClusterDepth.jl tutorial](https://www.s-ccs.de/ClusterDepth.jl/dev/tutorials/eeg/). -Let's create data from 20 subjects -````@example eeg +```@example eeg n_subjects = 20 data,events = UnfoldSim.predef_eeg(n_subjects;return_epoched=true) times = range(0,step=1/100,length=size(data,1)) -```` -Fit an UnfoldModel to each subject -````@example eeg +``` + + ```@raw html +
+ ``` + +## Fitting regression Models +Fit an UnfoldModel to each subject: +```@example eeg formula = @formula(0 ~ 1 + condition) models = map((d, ev) -> (fit(UnfoldModel, formula, DataFrame(ev), collect(d), times), ev.subject[1]), eachslice(data; dims=3), groupby(events, :subject)) -```` +nothing #hide + +``` now we can inspect the data easily, and extract the face-effect -````@example eeg +```@example eeg function add_subject!(df, s) df[!, :subject] .= s return df @@ -38,30 +57,30 @@ end allEffects = map((x) -> (effects(Dict(:condition => ["car", "face"]), x[1]), x[2]) |> (x) -> add_subject!(x[1], x[2]), models) |> e -> reduce(vcat, e) plot_erp(allEffects; mapping=(color=:condition, group=:subject)) -```` +``` extract the face-coefficient from the linear model -````@example eeg +```@example eeg allCoefs = map(m -> (coeftable(m[1]), m[2]) |> (x) -> add_subject!(x[1], x[2]), models) |> e -> reduce(vcat, e) plot_erp(allCoefs; mapping=(group=:subject, col=:coefname)) -```` +``` let's unstack the tidy-coef table into a matrix and put it to clusterdepth for clusterpermutation testing -````@example eeg +```@example eeg faceCoefs = allCoefs |> x -> subset(x, :coefname => x -> x .== "condition: face") erpMatrix = unstack(faceCoefs, :subject, :time, :estimate) |> x -> Matrix(x[:, 2:end])' |> collect summary(erpMatrix) -```` +``` ## Clusterdepth -````@example eeg +```@example eeg pvals = clusterdepth(erpMatrix; τ=quantile(TDist(n_subjects - 1), 0.95), nperm=5000); nothing #hide -```` +``` well - that was fast, less than a second for a cluster permutation test. not bad at all! @@ -70,18 +89,22 @@ Some plotting, and we add the identified cluster first calculate the ERP -````@example eeg +```@example eeg faceERP = groupby(faceCoefs, [:time, :coefname]) |> x -> combine(x, :estimate => mean => :estimate, :estimate => std => :stderror); nothing #hide -```` +``` + +put the significance into a dataframe-form -put the pvalues into a nicer format +```@example eeg +significant_regions_df = ClusterDepth.cluster(pvals .<= 0.05) |> x -> DataFrame(:from => times[x[1]], :to => times[x[1] .+ x[2]] , :coefname => "condition: face") +``` -````@example eeg -pvalDF = ClusterDepth.cluster(pvals .<= 0.05) |> x -> DataFrame(:from => x[1] ./ 250, :to => (x[1] .+ x[2]) ./ 250, :coefname => "condition: face") -plot_erp(faceERP; stderror=true, significance=pvalDF) -```` +and plot it! +```@example eeg +plot_erp(faceERP; stderror=true, significance=significant_regions_df) +``` Looks good to me! We identified the cluster :-) \ No newline at end of file From bdb39cc48306967360e9266f1a45e35afed37ac4 Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Tue, 18 Feb 2025 23:53:13 +0100 Subject: [PATCH 33/34] fix some irregularities --- docs/literate/tutorial/lmm_clusterdepth.jl | 15 ++++++++------- docs/src/tutorial/test_single_coefficient.md | 15 +++++++-------- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/docs/literate/tutorial/lmm_clusterdepth.jl b/docs/literate/tutorial/lmm_clusterdepth.jl index c60695d..f77333c 100644 --- a/docs/literate/tutorial/lmm_clusterdepth.jl +++ b/docs/literate/tutorial/lmm_clusterdepth.jl @@ -1,5 +1,5 @@ # !!! important -# This functionality is not tested. While it returns reasonable results, we haven't written any unit-tests, nor tested the type-1 error probability yet +# This functionality is not extensively tested. While it returns reasonable results, we haven't written any unit-tests, nor tested the type-1 error probability yet # ## 1. Simulate data @@ -57,10 +57,10 @@ data_e, events = UnfoldSim.simulate( UniformOnset(srate * 2, 10), PinkNoise(; noiselevel = 1); return_epoched = true, -) # 18 +) times = range(-0.1, 0.5, length = size(data_e, 1)) data_e = reshape(data_e, 1, size(data_e, 1), :) -#events.latency .+= repeat(range(0,length=size(data,2),step=size(data,1)),inner=size(events[events.subject.=="S01",:],1)) +nothing ##hide # ```@raw html # # ``` @@ -85,9 +85,9 @@ m = fit( # One solution are cluster permutation test, where we instead test for clustersizes of connected significant clusters. In "classical" two-stage testing, such a permutation test is straight forward. But for LMMs we have to think of something more clever, as it is not directly clear how to permute if both subject and item effects exist (you gonna break the relation between the two). We did that in `MixedModelsPermutations` and can apply this strategy to EEG data as well.` # select the fixed-effects coefficient to test (`stimtype`) -coefficient = 2 +coefficient = 2; -# call the permutatio test +# call the permutation test # !!! note # This interface is very likely to change in the future pvalue( @@ -95,6 +95,7 @@ pvalue( m, data_e, coefficient; - n_permutations = 10, + n_permutations = 20, clusterforming_threshold = 1.8, -) \ No newline at end of file +) + diff --git a/docs/src/tutorial/test_single_coefficient.md b/docs/src/tutorial/test_single_coefficient.md index 14e197a..2748aa9 100644 --- a/docs/src/tutorial/test_single_coefficient.md +++ b/docs/src/tutorial/test_single_coefficient.md @@ -8,10 +8,11 @@ We simulate a 1x2 design and perform cluster permutation testing via the `Cluste ## Simulate data Let's create data from 20 subjects - ```@raw html +```@raw html
Click to expand the simulation details - ``` +``` + ```@example eeg using Random using CairoMakie @@ -22,19 +23,17 @@ using Statistics using ClusterDepth using DataFrames using Distributions -``` - - -```@example eeg n_subjects = 20 data,events = UnfoldSim.predef_eeg(n_subjects;return_epoched=true) times = range(0,step=1/100,length=size(data,1)) + +nothing #hide ``` - ```@raw html +```@raw html
- ``` +``` ## Fitting regression Models Fit an UnfoldModel to each subject: From 1e4e12ff350c5beb8fb8665cd792a782fc432adf Mon Sep 17 00:00:00 2001 From: "behinger (s-ccs 001)" Date: Wed, 19 Feb 2025 09:37:34 +0100 Subject: [PATCH 34/34] minor fix --- docs/literate/tutorial/lmm_clusterdepth.jl | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/docs/literate/tutorial/lmm_clusterdepth.jl b/docs/literate/tutorial/lmm_clusterdepth.jl index f77333c..70e30ef 100644 --- a/docs/literate/tutorial/lmm_clusterdepth.jl +++ b/docs/literate/tutorial/lmm_clusterdepth.jl @@ -2,18 +2,22 @@ # This functionality is not extensively tested. While it returns reasonable results, we haven't written any unit-tests, nor tested the type-1 error probability yet +# ## 0. Setup +# If you want to do LMM + cluster permutations, you at least need these packages (more for the simulation below) +using UnfoldMixedModels +using MixedModelsPermutations, ClusterDepth +using UnfoldStats + # ## 1. Simulate data # This section can be skipped, if one already has (real) data that they want to analyse. + + # ```@raw html #
# Click to expand the simulation details # ``` using UnfoldSim -using UnfoldMixedModels - -using MixedModelsPermutations, ClusterDepth # both necessary to activate correct extension! -using UnfoldStats using StatsModels using Random