Repository: JuliaStats/MixedModels.jl (julia)

Issue:
Programmatically Creating RE Terms
This is at least partly a StatsModels.jl issue, partly a mixedmodels issue, and partly a julia/macro issue.

I think what you're trying to do is construct a formula from a list of variable names.  For more information on that, see the docs for statsmodels on constructing a formula at run-time: https://juliastats.org/StatsModels.jl/stable/formula/#Constructing-a-formula-programmatically-1.  You can do the fixed effects part easily:

```julia
julia> fix = term(1) + sum(term.("x$n" for n in 1:3))
1
x1(unknown)
x2(unknown)
x3(unknown)
```

The mixedmodels issue is that constructing a random-effects term at run-time using that technique is tricky, but it can be done.  The reason is that you have to construct the special `RandomEffectsTerm` used by mixedmodels, but that's assumed to happen after the schema has already been extracted (see this section of the statsmodels docs: https://juliastats.org/StatsModels.jl/stable/internals/#The-lifecycle-of-a-@formula-1, and this code from mixedmodels on how the random effects syntax with `|` is implemented: https://github.com/JuliaStats/MixedModels.jl/blob/master/src/randomeffectsterm.jl; note especially how the constructor for `RandomEffectsTerm` checks whether the grouping is an interaction or categorical term, which is only known after the schema is applied).  This is to my mind a legitimate quality of life issue.  At the moment you'll have to construct the concrete versions of all the terms and assemble the RandomEffectsTerm manually.  While this technically works (in that it will give you a formula you can fit with mixedmodels) it isn't going to give you any of the bumpers that mixedmodels normally provides (e.g., if you create multiple random effects terms this way, then you might end up with ... "interesting" (mis-specified) models), so it's definitely a "here be dragons" kinda thing that we should fix.  So with the caveat that I'm handing you a big foot-gun here, this is how I'd do it:

```julia
julia> ran = sum(term.("z$n" for n in 1:2))
z1(unknown)
z2(unknown)

julia> ran = apply_schema(sum(term.("z$n" for n in 1:2)), schema(dat), MixedModel)
z1(continuous)
z2(continuous)

julia> grouping = StatsModels.concrete_term(term(:g), dat, Dict(:g => Grouping()))
g(Grouping:5→5)

julia> reterm = MixedModels.RandomEffectsTerm(ran, grouping)
(z1 + z2 | g)

julia> formula = term(:y) ~ term(1) + fix + reterm
FormulaTerm
Response:
  y(unknown)
Predictors:
  1
  x1(unknown)
  x2(unknown)
  x3(unknown)
  (z1 + z2 | g)
```

This gets us most of the way there but it still doesn't work unfortunately.  The issue is that MixedModels tries to extract a schema and doesn't recognize that `RandomEffectsTerms` don't need them...we can duck patch that by defining `StatsModels.needs_schema(::RandomEffectsTerm) = false` but that's better fixed here or in statsmodels IMO.

In the interest of fixing this the right way, here are the missing pieces in MixedModels at the moment

- [ ] make `apply_schema(::FunctionTerm{typeof(|)}, ...)` create an "non-concrete" `RandomEffectsTerm`
- [ ] add additional `apply_schema(::RandomEffectsTerm, ...)` which makes it concrete, does the validity checks done in teh inner constructor right now
- [ ] add run-time methods for `Base.|(::AbstractTerm, ::AbstractTerm)` to create a `RandomEffectsTerm`, or if we don't want to do type piracy, use a different syntax (like `re` or `reterm` or something...)

Finally, the most superficial reason why your attempt doesn't work is a Julia/macro issue: `@formula(f)` only sees one thing as its argument, the Symbol `:f`.  Macros in julia can only do transformations of the surface syntax (symbols, exprs, strings, and literal numbers).  Reading the docs on macros in julia will be helpful, and people on slack/discourse/zulip can answer other questions in that direction.  If you're going to go this way (which is not recommended), then you need to use `eval` and `$f`.

_Originally posted by @kleinschmidt in https://github.com/JuliaStats/MixedModels.jl/issues/426#issuecomment-720564742_

Pull request description:
This fixes #462 and more or less follows the plan outlined there: it adds a
method for `Base.:|` for two terms which returns a `RandomEffectsTerm`, and adds
an additional layer of indirection for `apply_schema` to operate on the
"untyped" `RandomEffectsTerm`s created by that process.  I've also made the
analogous changes to `zerocorr` (although that runtime method was actually
already in place, it didn't, you know, work...)

To my mind, the one big wart here is that we're "claiming"
`Base.:|(::TermOrTerms, ::TermOrTerms)` which wouldn't play nicely with another
package that also wants to provide a different set of semantics for `|` (at
run-time).  I first thought that we could do the same trick that
FilePathsBase.jl does with `/` (creating a new `/` function that falls back to
`Base./` for arguments of types it doesn't own) but that didn't seem to work (I
got an "ignoring conflicting import" when doing `using MixedModels: |` with that
strategy).

Patch:
diff --git a/NEWS.md b/NEWS.md
index f5f391629..d9d7ede6a 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -9,6 +9,17 @@ MixedModels v4.0.0 Release Notes
   and make `MixedModelBootstrap` a subtype of it. Accordingly, rename the `bstr`
   field to `fits`. [#465]
 
+Run-time formula syntax
+-----------------------
+
+* It is now possible to construct `RandomEffectsTerm`s at run-time from `Term`s
+  (methods for `Base.|(::AbstractTerm, ::AbstractTerm)` added) [#470]
+* `RandomEffectsTerm`s can have left- and right-hand side terms that are
+  "non-concrete", and `apply_schema(::RandomEffectsTerm, ...)` works more like
+  other StatsModels.jl `AbstractTerm`s [#470]
+* Methods for `Base./(::AbstractTerm, ::AbstractTerm)` are added, allowing
+  nesting syntax to be used with `Term`s at run-time as well [#470]
+
 MixedModels v3.1.4 Release Notes
 ========================
 * [experimental] Additional convenience constructors for `LinearMixedModel` [#449]
diff --git a/docs/src/constructors.md b/docs/src/constructors.md
index f0aea5f92..b52403f03 100644
--- a/docs/src/constructors.md
+++ b/docs/src/constructors.md
@@ -27,13 +27,34 @@ dyestuff = DataFrame(MixedModels.dataset(:dyestuff))
 describe(dyestuff)
 ```
 
-### Models with simple, scalar random effects
+### The `@formula` language in Julia
 
-The formula language in *Julia* is similar to that in *R* which is based on (Wilkinson and Rogers [1973](https://dx.doi.org/10.2307/2346786)). In Julia a formula must be enclosed in a call to the `@formula` macro.
+MixedModels.jl builds on the the *Julia* formula language provided by [StatsModels.jl](https://juliastats.org/StatsModels.jl/stable/formula/), which is similar to the formula language in *R* and is also based on the notation from Wilkinson and Rogers ([1973](https://dx.doi.org/10.2307/2346786)). There are two ways to construct a formula in Julia.  The first way is to enclose the formula expression in the `@formula` macro:
 ```@docs
 @formula
 ```
 
+The second way is to combine `Term`s with operators like `+`, `&`, `~`, and others at "run time".  This is especially useful if you wish to create a formula from a list a variable names.  For instance, the following are equivalent:
+
+```@example Main
+@formula(y ~ 1 + a + b + a & b) == term(:y) ~ term(1) + term(:a) + term(:b) + term(:a) & term(:b)
+```
+
+MixedModels.jl provides additional formula syntax for representing *random-effects terms*.  Most importantly, `|` separates random effects and their grouping factors (as in the formula extension used by the *R* package [`lme4`](https://cran.r-project.org/web/packages/lme4/index.html).  Much like with the base formula language, `|` can be used within the `@formula` macro and to construct a formula programmatically:
+
+```@example Main
+@formula(y ~ 1 + a + b + (1 + a + b | g))
+```
+
+```@example Main
+terms = sum(term(t) for t in [1, :a, :b])
+group = term(:g)
+response = term(:y)
+response ~ terms + (terms | group)
+```
+
+### Models with simple, scalar random effects
+
 A basic model with simple, scalar random effects for the levels of `batch` (the batch of an intermediate product, in this case) is declared and fit as
 
 ```@example Main
@@ -535,4 +556,3 @@ fm4r = fit(MixedModel, @formula(diameter ~ 1+(1|plate)+(1|sample)),
 ```@example Main
 sum(leverage(fm4r))
 ```
-
diff --git a/src/randomeffectsterm.jl b/src/randomeffectsterm.jl
index f0b9ff350..e160ad467 100644
--- a/src/randomeffectsterm.jl
+++ b/src/randomeffectsterm.jl
@@ -3,43 +3,16 @@ abstract type AbstractReTerm <: AbstractTerm end
 struct RandomEffectsTerm <: AbstractReTerm
     lhs::StatsModels.TermOrTerms
     rhs::StatsModels.TermOrTerms
-    function RandomEffectsTerm(lhs, rhs)
-        if isempty(intersect(StatsModels.termvars(lhs), StatsModels.termvars(rhs)))
-            if !isa(
-                rhs,
-                Union{
-                    CategoricalTerm,
-                    InteractionTerm{<:NTuple{N,CategoricalTerm} where {N}},
-                },
-            )
-                throw(ArgumentError("blocking variables (those behind |) must be Categorical ($(rhs) is not)"))
-            end
-            new(lhs, rhs)
-        else
-            throw(ArgumentError("Same variable appears on both sides of |"))
-        end
-    end
 end
 
-function StatsModels.apply_schema(
-    t::FunctionTerm{typeof(/)},
-    sch::StatsModels.FullRank,
-    Mod::Type{<:MixedModel},
-)
-    if length(t.args_parsed) ≠ 2
-        throw(ArgumentError("malformed nesting term: $t (Exactly two arguments required"))
-    end
-
-    first, second = apply_schema.(t.args_parsed, Ref(sch), Mod)
-    
-    if !(typeof(first) <: CategoricalTerm)
-        throw(ArgumentError("nesting terms requires categorical grouping term, got $first.  Manually specify $first as `CategoricalTerm` in hints/contrasts"))
-    end
-
-    return first + fulldummy(first) & second
-end
+# TODO: consider overwriting | with our own function that can be
+# imported with (a la FilePathsBase.:/)
+# using MixedModels: |
+# to avoid conflicts with definitions in other packages...
+Base.:|(a::StatsModels.TermOrTerms, b::StatsModels.TermOrTerms) = RandomEffectsTerm(a, b)
 
-RandomEffectsTerm(lhs, rhs::NTuple{2,AbstractTerm}) =
+# expand (lhs | a + b) to (lhs | a) + (lhs | b)
+RandomEffectsTerm(lhs::StatsModels.TermOrTerms, rhs::NTuple{2,AbstractTerm}) =
     (RandomEffectsTerm(lhs, rhs[1]), RandomEffectsTerm(lhs, rhs[2]))
 
 Base.show(io::IO, t::RandomEffectsTerm) = print(io, "($(t.lhs) | $(t.rhs))")
@@ -49,6 +22,9 @@ function StatsModels.termvars(t::RandomEffectsTerm)
     vcat(StatsModels.termvars(t.lhs), StatsModels.termvars(t.rhs))
 end
 
+StatsModels.terms(t::RandomEffectsTerm) = union(StatsModels.terms(t.lhs), StatsModels.terms(t.rhs))
+
+# | in MixedModel formula -> RandomEffectsTerm
 function StatsModels.apply_schema(
     t::FunctionTerm{typeof(|)},
     schema::MultiSchema{StatsModels.FullRank},
@@ -56,8 +32,30 @@ function StatsModels.apply_schema(
 )
     lhs, rhs = t.args_parsed
 
+    isempty(intersect(StatsModels.termvars(lhs), StatsModels.termvars(rhs))) ||
+        throw(ArgumentError("Same variable appears on both sides of |"))
+
+    apply_schema(RandomEffectsTerm(lhs, rhs), schema, Mod)
+end
+
+# allowed types (or tuple thereof) for blocking variables (RHS of |):
+const GROUPING_TYPE = Union{<:CategoricalTerm, <:InteractionTerm{<:NTuple{N,CategoricalTerm} where {N}}}
+check_re_group_type(term::GROUPING_TYPE) = true
+check_re_group_type(terms::Tuple{Vararg{<:GROUPING_TYPE}}) = true
+check_re_group_type(x) = false
+
+# make a potentially untyped RandomEffectsTerm concrete
+function StatsModels.apply_schema(
+    t::RandomEffectsTerm,
+    schema::MultiSchema{StatsModels.FullRank},
+    Mod::Type{<:MixedModel}
+)
+    lhs, rhs = t.lhs, t.rhs
+
+    # get a schema that's specific for the grouping (RHS), creating one if needed
     schema = get!(schema.subs, rhs, StatsModels.FullRank(schema.base.schema))
 
+    # handle intercept in LHS (including checking schema for intercept in another term)
     if (
         !StatsModels.hasintercept(lhs) &&
         !StatsModels.omitsintercept(lhs) &&
@@ -66,7 +64,13 @@ function StatsModels.apply_schema(
     )
         lhs = InterceptTerm{true}() + lhs
     end
+
     lhs, rhs = apply_schema.((lhs, rhs), Ref(schema), Mod)
+
+    # check whether grouping terms are categorical or interaction of categorical
+    check_re_group_type(rhs) ||
+        throw(ArgumentError("blocking variables (those behind |) must be Categorical ($(rhs) is not)"))
+    
     RandomEffectsTerm(MatrixTerm(lhs), rhs)
 end
 
@@ -98,7 +102,6 @@ function StatsModels.modelcols(t::RandomEffectsTerm, d::NamedTuple)
     )
 end
 
-
 # extract vector of refs from ranef grouping term and data
 function _ranef_refs(grp::CategoricalTerm, d::NamedTuple)
     invindex = grp.contrasts.invindex
@@ -117,6 +120,25 @@ function _ranef_refs(
     refs, uniques
 end
 
+# TODO: split this off into a RegressionFormula packge?
+Base.:/(a::AbstractTerm, b::AbstractTerm) = a + a & b
+function StatsModels.apply_schema(
+    t::FunctionTerm{typeof(/)},
+    sch::StatsModels.FullRank,
+    Mod::Type{<:MixedModel},
+)
+    if length(t.args_parsed) ≠ 2
+        throw(ArgumentError("malformed nesting term: $t (Exactly two arguments required"))
+    end
+
+    first, second = apply_schema.(t.args_parsed, Ref(sch), Mod)
+
+    if !(typeof(first) <: CategoricalTerm)
+        throw(ArgumentError("nesting terms requires categorical grouping term, got $first.  Manually specify $first as `CategoricalTerm` in hints/contrasts"))
+    end
+
+    return first + fulldummy(first) & second
+end
 
 # add some syntax to manually promote to full dummy coding
 fulldummy(t::AbstractTerm) =
@@ -169,6 +191,10 @@ Remove correlations between random effects in `term`.
 """
 zerocorr(x) = ZeroCorr(x)
 
+# for schema extraction (from runtime-created zerocorr)
+StatsModels.terms(t::ZeroCorr) = StatsModels.terms(t.term)
+StatsModels.termvars(t::ZeroCorr) = StatsModels.termvars(t.term)
+
 function StatsModels.apply_schema(
     t::FunctionTerm{typeof(zerocorr)},
     sch::MultiSchema,
@@ -177,4 +203,12 @@ function StatsModels.apply_schema(
     ZeroCorr(apply_schema(t.args_parsed..., sch, Mod))
 end
 
+function StatsModels.apply_schema(
+    t::ZeroCorr,
+    sch::MultiSchema,
+    Mod::Type{<:MixedModel},
+)
+    ZeroCorr(apply_schema(t.term, sch, Mod))
+end
+
 StatsModels.modelcols(t::ZeroCorr, d::NamedTuple) = zerocorr!(modelcols(t.term, d))


Test patch:
diff --git a/test/FactorReTerm.jl b/test/FactorReTerm.jl
index 7dadd0335..68e83e0bf 100644
--- a/test/FactorReTerm.jl
+++ b/test/FactorReTerm.jl
@@ -99,6 +99,55 @@ end
         form = apply_schema(f, schema(f, slp, contrasts), LMM)
         @test StatsModels.termvars(form.rhs) == [:days, :subj]
     end
+
+    @testset "Runtime construction of random effects terms" begin
+        # operator precedence and basic terms:
+        @test term(:a) | term(:b) isa RandomEffectsTerm
+        @test term(1) + term(:a) | term(:b) isa RandomEffectsTerm
+        @test term(1) + term(:a) + term(:a) & term(:c) | term(:b) isa RandomEffectsTerm
+        
+        # sleep study data:
+        r, d, s, one = term.((:reaction, :days, :subj, 1))
+
+        f1 = @formula(reaction ~ 1 + (1 + days | subj))
+        f2 = r ~ one + (one + d | s)
+        @test f2.rhs[end] isa RandomEffectsTerm
+        ff1 = apply_schema(f1, schema(slp), LMM)
+        ff2 = apply_schema(f2, schema(slp), LMM)
+        # equality of RE terms not defined so check that they generate same modelcols
+        @test modelcols(ff1.rhs[end], slp) == modelcols(ff2.rhs[end], slp)
+        
+        m1 = fit(LMM, f1, slp)
+        m2 = fit(LMM, f2, slp)
+        @test all(m1.λ .== m2.λ)
+
+        @test StatsModels.terms(f2.rhs[end]) == [one, d, s]
+        @test StatsModels.termvars(f2.rhs[end]) == [d.sym, s.sym]
+    end
+
+    @testset "Runtime construction of ZeroCorr" begin
+        r, d, s, one = term.((:reaction, :days, :subj, 1))
+
+        f1 = @formula(reaction ~ 1 + zerocorr(1 + days | subj))
+        f2 = r ~ one + zerocorr(one + d | s)
+        @test f2.rhs[end] isa MixedModels.ZeroCorr
+        ff1 = apply_schema(f1, schema(slp), LMM)
+        ff2 = apply_schema(f2, schema(slp), LMM)
+        # equality of RE terms not defined so check that they generate same modelcols
+        mc1 = modelcols(ff1.rhs[end], slp)
+        mc2 = modelcols(ff2.rhs[end], slp)
+
+        # test that zerocorr actually worked
+        @test mc1.inds == mc2.inds == [1, 4]
+        
+        m1 = fit(LMM, f1, slp)
+        m2 = fit(LMM, f2, slp)
+        @test all(m1.λ .== m2.λ)
+
+        @test StatsModels.terms(f2.rhs[end]) == [one, d, s]
+        @test StatsModels.termvars(f2.rhs[end]) == [d.sym, s.sym]
+    end
+    
 end
 
 @testset "Categorical Blocking Variable" begin
@@ -181,10 +230,31 @@ end
         @test modelcols(f2.rhs, d2) == [ones(20) d2.b .== :Y (d2.b .== :X).*d2.a (d2.b .== :Y).*d2.a]
         @test coefnames(f2.rhs) == ["(Intercept)", "b: Y", "b: X & a", "b: Y & a"]
 
+        # runtime
+        fr2 = term(1) + term(:b) / term(:a)
+        @test length(fr2) == 3
+        @test fr2[1] isa ConstantTerm
+        @test fr2[2] isa Term
+        @test fr2[3] isa InteractionTerm
+        frf2 = apply_schema(term(0) ~ fr2, schema(d2), MixedModel)
+        @test modelcols(frf2.rhs, d2) == [ones(20) d2.b .== :Y (d2.b .== :X).*d2.a (d2.b .== :Y).*d2.a]
+        @test coefnames(frf2.rhs) == ["(Intercept)", "b: Y", "b: X & a", "b: Y & a"]
+
+        # check promotion
         f3 = apply_schema(@formula(0 ~ 0 + b/a), schema(d2), MixedModel)
         @test modelcols(f3.rhs, d2) == [d2.b .== :X d2.b .== :Y (d2.b .== :X).*d2.a (d2.b .== :Y).*d2.a]
         @test coefnames(f3.rhs) == ["b: X", "b: Y", "b: X & a", "b: Y & a"]
 
+        # runtime:
+        fr3 = term(0) + term(:b) / term(:a)
+        @test length(fr3) == 3
+        @test fr3[1] isa ConstantTerm
+        @test fr3[2] isa Term
+        @test fr3[3] isa InteractionTerm
+        ffr3 = apply_schema(term(0) ~ fr3, schema(d2), MixedModel)
+        @test modelcols(ffr3.rhs, d2) == [d2.b .== :X d2.b .== :Y (d2.b .== :X).*d2.a (d2.b .== :Y).*d2.a]
+        @test coefnames(ffr3.rhs) == ["b: X", "b: Y", "b: X & a", "b: Y & a"]
+
         # errors for continuous grouping
         @test_throws ArgumentError apply_schema(@formula(0 ~ 1 + a/b), schema(d2), MixedModel)
 
@@ -195,10 +265,12 @@ end
         psts = dataset("pastes")
         m = fit(MixedModel, @formula(strength ~ 1 + (1|batch/cask)), psts)
         m2 = fit(MixedModel, @formula(strength ~ 1 + (1|batch) + (1|batch&cask)), psts)
+        mr = fit(MixedModel, term(:strength) ~ term(1) + (term(1)|term(:batch)/term(:cask)), psts)
+        m2r = fit(MixedModel, term(:strength) ~ term(1) + (term(1)|term(:batch)) + (term(1)|term(:batch)&term(:cask)), psts)
 
-        @test fnames(m) == fnames(m2) == (Symbol("batch & cask"), :batch)
-        @test m.λ == m2.λ
-        @test deviance(m) == deviance(m2)
+        @test fnames(m) == fnames(m2) == fnames(mr) == fnames(m2r) == (Symbol("batch & cask"), :batch)
+        @test m.λ == m2.λ == mr.λ == m2r.λ
+        @test deviance(m) == deviance(m2) == deviance(mr) == deviance(m2r)
     end
 
     @testset "multiple terms with same grouping" begin