Repository: JuliaStats/MixedModels.jl (julia)

Issue:
Automatic conversion of grouping variables to categorical for @formula but not term() syntax
MWE
```julia
using MixedModels, StatsModels
df = (x=randn(1000), y=randn(1000), g=rand([0,1,2,3],1000))
m = fit(LinearMixedModel, @formula(y ~ x + (1|g)), df) # no error
m1 = fit(LinearMixedModel, term(:y) ~ term(:x) + (term(1)|term(:g)), df)
# ERROR: ArgumentError: blocking variables (those behind |) must be Categorical (g is not)
```

Maybe this just needs to be a change in documentation, but it's confusing that the `@formula` syntax does conversion but the `term(..)` syntax errors.

Pull request:
- [x] add entry in NEWS.md
- [x] after opening this PR, add a reference and run `docs/NEWS-update.jl` to update the cross-references.
- [x] I've bumped the version appropriately

closes #836

Patch:
diff --git a/.github/workflows/minimum.yml b/.github/workflows/minimum.yml
index 4e9460969..a4886b174 100644
--- a/.github/workflows/minimum.yml
+++ b/.github/workflows/minimum.yml
@@ -21,7 +21,7 @@ jobs:
     strategy:
       matrix:
         julia-version: [min]
-        os: [ubuntu-22.04, macos-14, windows-2019]
+        os: [ubuntu-22.04, macos-14, windows-2022]
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
diff --git a/NEWS.md b/NEWS.md
index c45e8b81f..bdcfb58a3 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,7 @@
+MixedModels v4.36.0 Release Notes
+==============================
+- Automatic application of grouping contrasts now works for `RandomEffectsTerm`s constructed programmatically, e.g. with `(term(1) | term(:g))`. [#836]
+
 MixedModels v4.35.2 Release Notes
 ==============================
 - Update to the `show(::IO, ::MIME, x)` methods for consistent ordering across Julia versions of variance components without associated fixed effects. This may result in change of ordering on existing Julia versions, but ordering will now match the ordering in the default REPL display (and reflect the model internal ordering), [#829]
@@ -642,3 +646,4 @@ Package dependencies
 [#825]: https://github.com/JuliaStats/MixedModels.jl/issues/825
 [#828]: https://github.com/JuliaStats/MixedModels.jl/issues/828
 [#829]: https://github.com/JuliaStats/MixedModels.jl/issues/829
+[#836]: https://github.com/JuliaStats/MixedModels.jl/issues/836
diff --git a/Project.toml b/Project.toml
index 715fc587a..d7b2a2ffe 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "MixedModels"
 uuid = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 author = ["Phillip Alday <me@phillipalday.com>", "Douglas Bates <dmbates@gmail.com>"]
-version = "4.35.2"
+version = "4.36.0"
 
 [deps]
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"
@@ -57,6 +57,7 @@ PooledArrays = "0.5, 1"
 PrecompileTools = "1"
 ProgressMeter = "1.7"
 Random = "1"
+RegressionFormulae = "0.1.3"
 SparseArrays = "1"
 StableRNGs = "0.1, 1"
 StaticArrays = "0.11, 0.12, 1"
@@ -78,9 +79,10 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 PRIMA = "0a7d04aa-8ac2-47b3-b7a7-9dbd6ad661ed"
+RegressionFormulae = "545c379f-4ec2-4339-9aea-38f2fb6a8ba2"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "DataFrames", "ExplicitImports", "InteractiveUtils", "PRIMA", "StableRNGs", "Suppressor", "Test"]
+test = ["Aqua", "DataFrames", "ExplicitImports", "InteractiveUtils", "PRIMA", "RegressionFormulae", "StableRNGs", "Suppressor", "Test"]
diff --git a/README.md b/README.md
index 682501384..66e163607 100644
--- a/README.md
+++ b/README.md
@@ -56,7 +56,7 @@ Typical distribution forms are _Bernoulli_ for binary data or _Poisson_ for coun
 |Linux   | Ubuntu 22.04  | x64    |current release |
 |Linux   | Ubuntu 22.04  | x64    |nightly         |
 |macOS   | Sonoma 14     | aarm64 |v1.10           |
-|Windows | Server 2019   | x64    |v1.10           |
+|Windows | Server 2022   | x64    |v1.10           |
 
 Note that previous releases still support older Julia versions.
 
diff --git a/src/MixedModels.jl b/src/MixedModels.jl
index 561e970bd..b991b4ac8 100644
--- a/src/MixedModels.jl
+++ b/src/MixedModels.jl
@@ -40,7 +40,7 @@ using StatsFuns: log2π, normccdf
 using StatsModels: StatsModels, AbstractContrasts, AbstractTerm, CategoricalTerm
 using StatsModels: ConstantTerm, DummyCoding, EffectsCoding, FormulaTerm, FunctionTerm
 using StatsModels: HelmertCoding, HypothesisCoding, InteractionTerm, InterceptTerm
-using StatsModels: MatrixTerm, SeqDiffCoding, TableRegressionModel, Term
+using StatsModels: MatrixTerm, SeqDiffCoding, TableRegressionModel
 using StatsModels: apply_schema, drop_term, formula, lrtest, modelcols, @formula
 using StructTypes: StructTypes
 using Tables: Tables, columntable
diff --git a/src/randomeffectsterm.jl b/src/randomeffectsterm.jl
index f5ef4a6dd..c7fdacdf0 100644
--- a/src/randomeffectsterm.jl
+++ b/src/randomeffectsterm.jl
@@ -33,15 +33,20 @@ end
 
 schema(t, data, hints) = StatsModels.schema(t, data, hints)
 
-function schema(t::FunctionTerm{typeof(|)}, data, hints::Dict{Symbol})
-    sch = schema(t.args[1], data, hints)
-    vars = StatsModels.termvars.(t.args[2])
+function schema(t::AbstractReTerm, data, hints::Dict{Symbol})
+    sch = schema(t.lhs, data, hints)
+    vars = StatsModels.termvars.(t.rhs)
     # in the event that someone has x|x, then the Grouping()
-    # gets overwrriten by the broader schema BUT
+    # gets overwritten by the broader schema BUT
     # that doesn't matter because we detect and throw an error
     # for that in apply_schema
     grp_hints = Dict(rr => Grouping() for rr in vars)
-    return merge(schema(t.args[2], data, grp_hints), sch)
+    return merge(schema(t.rhs, data, grp_hints), sch)
+end
+
+function schema(t::FunctionTerm{typeof(|)}, data, hints::Dict{Symbol})
+    re = RandomEffectsTerm(t.args[1], t.args[2])
+    return schema(re, data, hints)
 end
 
 is_randomeffectsterm(::Any) = false
@@ -101,7 +106,8 @@ function StatsModels.apply_schema(
     )
         lhs = InterceptTerm{true}() + lhs
     end
-    lhs, rhs = apply_schema.((lhs, rhs), Ref(schema), Mod)
+    lhs = apply_schema(lhs, schema, Mod)
+    rhs = apply_schema(rhs, schema, Mod)
 
     # check whether grouping terms are categorical or interaction of categorical
     check_re_group_type(rhs) || throw(
diff --git a/src/schema.jl b/src/schema.jl
index aae01520b..e08ac424b 100644
--- a/src/schema.jl
+++ b/src/schema.jl
@@ -56,7 +56,7 @@ end
 Find and apply the schema for f in a way that automatically uses `Grouping()`
 contrasts when appropriate.
 
-!!! warn
+!!! warning
     This is an internal method.
 """
 function schematize(f, tbl, contrasts::Dict{Symbol}, Mod=LinearMixedModel)
@@ -64,11 +64,21 @@ function schematize(f, tbl, contrasts::Dict{Symbol}, Mod=LinearMixedModel)
     # also we want this to be a vector so we can sort later
     rhs = f.rhs isa AbstractTerm ? [f.rhs] : collect(f.rhs)
     fe = filter(!is_randomeffectsterm, rhs)
+    # flatten any MatrixTerms
+    fe = mapreduce(StatsModels.terms, vcat, fe; init=AbstractTerm[])
     # init with lhs so we don't need an extra merge later
     # and so that things work even when we have empty fixed effects
     init = schema(f.lhs, tbl, contrasts)
     sch_fe = mapfoldl(merge, fe; init) do tt
-        return schema(tt, tbl, contrasts)
+        sch = schema(tt, tbl, contrasts)
+        # this happens when the term is already concrete
+        # we want to force the schema to populated so that it's not overridden
+        # by any associated RE below
+        if isempty(sch.schema)
+            sch = StatsModels.Schema(StatsModels.term(string(tt)) => tt)
+        end
+
+        return sch
     end
     re = filter(is_randomeffectsterm, rhs)
     sch_re = mapfoldl(merge, re; init) do tt
diff --git a/src/utilities.jl b/src/utilities.jl
index 413413b7c..0fa46f43a 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -13,13 +13,10 @@ function _abstractify_grouping(f::FormulaTerm)
     re = map(re) do trm
         if trm.rhs isa InteractionTerm
             rhs = mapreduce(&, trm.rhs.terms) do tt
-                # how to define Grouping() for interactions on the RHS?
-                # contr[tt.sym] = Grouping()
-                return Term(tt.sym)
+                return StatsModels.term(tt.sym)
             end
         else
-            contr[trm.rhs.sym] = Grouping()
-            rhs = Term(trm.rhs.sym)
+            rhs = StatsModels.term(trm.rhs.sym)
         end
         return trm.lhs | rhs
     end


Test patch:
diff --git a/test/grouping.jl b/test/grouping.jl
index 834f45282..b7fd91985 100644
--- a/test/grouping.jl
+++ b/test/grouping.jl
@@ -34,11 +34,14 @@ end
     d = (; y=rand(100),
         x=rand('A':'Z', 100),
         z=rand('A':'Z', 100),
-        grp=rand('a':'z', 100))
+        grp=rand(1:26, 100)) # we want this to be numeric so that we don't get past categorical checks by default
     contrasts = Dict{Symbol, Any}()
 
     @testset "blocking variables are grouping" for f in [@formula(y ~ 1 + x + (1|grp)),
-                                                         @formula(y ~ 1 + x + zerocorr(1|grp))]
+                                                         @formula(y ~ 1 + x + zerocorr(1|grp)),
+                                                         term(:y) ~ term(:x) + (term(1)|(term(:grp))),
+                                                         term(:y) ~ term(:x) + zerocorr(term(1)|(term(:grp)))
+                                                        ]
         fsch = schematize(f, d, contrasts)
         fe = fsch.rhs[1]
         x = last(fe.terms)
@@ -49,7 +52,10 @@ end
     end
 
     @testset "FE contrasts take priority" for f in [@formula(y ~ 1 + x + (1|x)),
-                                                    @formula(y ~ 1 + x + zerocorr(1|x))]
+                                                    @formula(y ~ 1 + x + zerocorr(1|x)),
+                                                    term(:y) ~ term(:x) + (term(1)|(term(:x))),
+                                                    term(:y) ~ term(:x) + zerocorr(term(1)|(term(:x)))
+                                                    ]
         fsch = schematize(f, d, contrasts)
         fe = fsch.rhs[1]
         x = last(fe.terms)
@@ -67,8 +73,11 @@ end
         @test grp.contrasts isa ContrastsMatrix{EffectsCoding}
     end
 
-    @testset "Nesting and interactions" for f in [@formula(y ~ 1 + x + (1|grp/z))]
+    @testset "Nesting and interactions" for f in [@formula(y ~ 1 + x + (1 | grp/z))]
         # XXX zerocorr(1|grp/z) doesn't work!
+        # XXX programmatic form doesn't work: term(:y) ~ term(:x) + (1 | term(:grp) / term(:z))
+        #  - we don't define ~pirate~ define the relevant methods here
+        #  - and RegressionFormulae.jl errors here as well
         fsch = schematize(f, d, contrasts)
         fe = fsch.rhs[1]
         x = last(fe.terms)
diff --git a/test/runtests.jl b/test/runtests.jl
index 307e9aec2..98bb8c06c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,6 +2,10 @@ using Aqua
 using ExplicitImports
 using GLM
 using MixedModels
+# we don't actually use RegressionFormulae
+# but we want to make sure we're not at odds with each other
+# in our term privateering
+using RegressionFormulae
 using Test
 
 import InteractiveUtils: versioninfo