Merge pull request #14 from unfoldtoolbox/bug-#13_1x3

Fix 1x3 categorical condensing

Author: behinger

Project.toml

@@ -1,7 +1,7 @@
 name = "UnfoldMixedModels"
 uuid = "019ae9e0-8363-565c-86e5-97a5a2fe84f4"
 authors = ["Benedikt Ehinger <benedikt.ehinger@vis.uni-stuttgart.de>"]
-version = "0.1.1"
+version = "0.1.2"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"

src/condense.jl

@@ -10,15 +10,26 @@ end
 MixedModels.tidyβ(m::Union{UnfoldLinearMixedModel,UnfoldLinearMixedModelContinuousTime}) =
     MixedModels.tidyβ(modelfit(m))
 
+
+"""
+helper function because `coefnames` returns an array only if number of coefs is larger than 1
+"""
+function _coefnames(t)
+    c = coefnames(t)
+    return isa(c, Vector) ? c : [c]
+end
+
+
+
 """
     random_effect_groupings(t::MixedModels.AbstractReTerm)
 Returns the random effect grouping term (rhs), similar to coefnames, which returns the left hand sides
 """
-random_effect_groupings(t::AbstractTerm) = repeat([nothing], length(t.terms))
+random_effect_groupings(t::AbstractTerm) = repeat([nothing], length(_coefnames(t.terms)))
 random_effect_groupings(t::Unfold.TimeExpandedTerm) =
     repeat(random_effect_groupings(t.term), length(Unfold.colnames(t.basisfunction)))
 random_effect_groupings(t::MixedModels.AbstractReTerm) =
-    repeat([t.rhs.sym], length(t.lhs.terms))
+    repeat([t.rhs.sym], length(_coefnames(t.lhs)))
 
 random_effect_groupings(f::FormulaTerm) = vcat(random_effect_groupings.(f.rhs)...)
 random_effect_groupings(t::Vector) = vcat(random_effect_groupings.(t)...)
@@ -32,16 +43,18 @@ function reorder_tidyσs(t, f)
     #@debug typeof(f)
     # get the order from the formula, this is the target
     f_order = random_effect_groupings(f) # formula order
-    #@debug f_order
+    @debug f_order
     f_order = vcat(f_order...)
-    #@debug f_order
+    @debug f_order
 
     # find the fixefs
     fixef_ix = isnothing.(f_order)
 
 
 
     f_order = string.(f_order[.!fixef_ix])
+    @debug fixef_ix
+    @debug coefnames(f)
 
     f_name = vcat(coefnames(f)...)[.!fixef_ix]
 
@@ -50,8 +63,8 @@ function reorder_tidyσs(t, f)
     t_name = [string(i.column) for i in t if i.iter == 1]
 
     # combine for formula and tidy output the group + the coefname
-    #@debug f_order
-    #@debug f_name
+    @debug "f" f_order f_name
+    @debug "t" t_order t_name
     f_comb = f_order .* f_name
     t_comb = t_order .* t_name
 
@@ -68,7 +81,7 @@ function reorder_tidyσs(t, f)
     # repeat and build the index for all timepoints
     reorder_ix_all = repeat(reorder_ix, length(t) ÷ length(reorder_ix))
     for k = 1:length(reorder_ix):length(t)
-        reorder_ix_all[k:k+length(reorder_ix)-1] .+= (k - 1)
+        reorder_ix_all[k:(k+length(reorder_ix)-1)] .+= (k - 1)
     end
 
     return t[reorder_ix_all]

test/test-fit.jl

@@ -24,7 +24,7 @@ end
     subj_idx = [parse(Int, split(string(s), 'S')[2]) for s in evts.subject]
     evts.latency .+= size(data, 1) .* (subj_idx .- 1)
 
-
+    evts.C = rand(StableRNG(1), ["a", "b", "c"], size(evts, 1))
     data = reshape(data, 1, :)
     #append!(data, zeros(1000))
     data = vcat(data, data)
@@ -62,7 +62,10 @@ end
     )
     df = Unfold.coeftable(m_mum)
     @test isapprox(
-        df[(df.channel.==1).&(df.coefname.=="A: a_small").&(df.time.==0.0), :estimate],
+        df[
+            (df.channel .== 1) .& (df.coefname .== "A: a_small") .& (df.time .== 0.0),
+            :estimate,
+        ],
         [-0.02, 0.054],
         atol = 0.01,
     )
@@ -81,7 +84,10 @@ end
     )
     df = coeftable(m_mum)
     @test isapprox(
-        df[(df.channel.==1).&(df.coefname.=="A: a_small").&(df.time.==0.0), :estimate],
+        df[
+            (df.channel .== 1) .& (df.coefname .== "A: a_small") .& (df.time .== 0.0),
+            :estimate,
+        ],
         [0.031, 0.05],
         atol = 0.1,
     )
@@ -101,7 +107,10 @@ end
     )
     df = coeftable(m_tum)
     @test isapprox(
-        df[(df.channel.==1).&(df.coefname.=="A: a_small").&(df.time.==0.0), :estimate],
+        df[
+            (df.channel .== 1) .& (df.coefname .== "A: a_small") .& (df.time .== 0.0),
+            :estimate,
+        ],
         [-0.03, 0.064],
         atol = 0.1,
     )
@@ -120,13 +129,12 @@ end
 
 
     evts.subjectB = evts.subject
-    evts1 = evts[evts.A.=="a_small", :]
-    evts2 = evts[evts.A.=="a_big", :]
+    evts1 = evts[evts.A .== "a_small", :]
+    evts2 = evts[evts.A .== "a_big", :]
 
-    f0_lmm = @formula 0 ~ 1 + B + (1 | subject) + (1 | subjectB)
-    @time m_tum = coeftable(
-        fit(UnfoldModel, f0_lmm, evts, data, basisfunction; show_progress = false),
-    )
+    f0_lmm = @formula 0 ~ 1 + B + (1 | subject) + (1 + C | subjectB)
+    @time m = fit(UnfoldModel, f0_lmm, evts, data, basisfunction; show_progress = false)
+    @time m_tum = coeftable(m)
 
 
     f1_lmm = @formula 0 ~ 1 + B + (1 | subject)
@@ -147,7 +155,10 @@ end
     df = coeftable(r)
 
     @test isapprox(
-        df[(df.channel.==1).&(df.coefname.=="B: b_tiny").&(df.time.==0.0), :estimate],
+        df[
+            (df.channel .== 1) .& (df.coefname .== "B: b_tiny") .& (df.time .== 0.0),
+            :estimate,
+        ],
         [0.65, 0.69],
         rtol = 0.1,
     )
@@ -163,6 +174,35 @@ end
         ),
     )
 
+
+    #----
+    # # test #13, 2x3 design
+
+    f = @formula 0 ~ 1 + A + C + (1 + A + C | subject)
+    #f  = @formula 0~1 + (1|subject)
+
+
+
+    # cut the data into epochs
+    # TODO This ignores subject bounds
+    data_e, times = Unfold.epoch(data = data, tbl = evts, τ = (-1.0, 1.9), sfreq = 10)
+    data_missing_e, times =
+        Unfold.epoch(data = data_missing, tbl = evts, τ = (-1.0, 1.9), sfreq = 10)
+    evts_e, data_e = Unfold.drop_missing_epochs(copy(evts), data_e)
+    evts_missing_e, data_missing_e = Unfold.drop_missing_epochs(copy(evts), data_missing_e)
+
+    ######################
+    ##  Mass Univariate Mixed
+    @time m_mum = fit(
+        UnfoldModel,
+        f,
+        evts_e,
+        data_e,
+        times,
+        contrasts = Dict(:A => EffectsCoding(), :B => EffectsCoding()),
+        show_progress = false,
+    )
+    df = Unfold.coeftable(m_mum)
 end
 ## Condense check for multi channel, multi
 @testset "LMM multi channel, multi basisfunction" begin
@@ -198,6 +238,22 @@ end
 
     @test all(last(.!isnothing.(res.group), 8))
     @test all(last(res.coefname, 8) .== "(Intercept)")
+
+    # test more complex formulas
+    fA0 = @formula (0 ~ 1 + zerocorr(1 + C | subject))
+    fA1 = @formula (0 ~ 1 + B + C + zerocorr(1 + C | subject2))
+    evts.C = rand(StableRNG(1), ["a", "b", "c"], size(evts, 1))
+    m = fit(
+        UnfoldModel,
+        ["a_small" => (fA0, bA0), "a_big" => (fA1, bA1)],
+        evts,
+        data;
+        eventcolumn = "A",
+        show_progress = false,
+    )
+
+    res = coeftable(m)
+
 end
 
 
@@ -216,17 +272,15 @@ end
         [
             Any => (
                 @formula(
-                    0 ~
-                        1 + A + B + zerocorr(1 + B + A | subject) + zerocorr(1 + B | item)
+                    0 ~ 1 + A + B + zerocorr(1 + B + A | subject) + zerocorr(1 + B | item)
                 ),
                 range(0, 1, length = size(data, 1)),
             ),
         ],
         [
             Any => (
                 @formula(
-                    0 ~
-                        1 + A + B + zerocorr(1 + A + B | subject) + zerocorr(1 + B | item)
+                    0 ~ 1 + A + B + zerocorr(1 + A + B | subject) + zerocorr(1 + B | item)
                 ),
                 range(0, 1, length = size(data, 1)),
             ),