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 1/8] 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 2/8] 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 3/8] 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 4/8] 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 5/8] 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 6/8] 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 7/8] 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 8/8] 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