# Should model matrix for nested factors be full-rank?

**URL:** <https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733>\
**Category:** Statistics\
**Created:** [July 25, 2022, 7:17am UTC](https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733 "2022-07-25T07:17:03Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [July 25, 2022, 7:17am UTC](https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733/1 "2022-07-25T07:17:03Z")

</div>

In some cases when levels of one factor all “inside” levels of another we have not full-rank model matrix from StatsModels. For example - calssical bioequivlence crossover study (2X2) with factors: `Var ~ Formulation + Sequence + Period + Subject` - levels of `Subject` distributed in `Sequnce`. So we get result bellow.

Should we get operator for `@formula` to make coding `inside` each level of top level factor?

```julia
:(log(Cmax)) ~ 1 + Sequence + Period + Formulation + Subject

Coefficients:
────────────────────────────────────────────────────────────────────────────────────
                     Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95%
────────────────────────────────────────────────────────────────────────────────────
(Intercept) 6.46193 0.077827 83.03 <1e-28 6.30052 6.62333
Sequence: tr -0.0143149 0.105746 -0.14 0.8935 -0.233619 0.204989
Period: 2 0.039951 0.0305263 1.31 0.2041 -0.0233566 0.103259
Formulation: t 0.0358771 0.0305263 1.18 0.2524 -0.0274304 0.0991847
Subject: 2 0.518319 0.105746 4.90 <1e-04 0.299015 0.737623
Subject: 3 0.301846 0.105746 2.85 0.0092 0.082542 0.52115
Subject: 4 0.0409981 0.105746 0.39 0.7020 -0.178306 0.260302
Subject: 5 -0.596241 0.105746 -5.64 <1e-04 -0.815545 -0.376937
Subject: 6 -0.950327 0.105746 -8.99 <1e-08 -1.16963 -0.731023
Subject: 7 0.115056 0.105746 1.09 0.2884 -0.104247 0.33436
Subject: 8 0.61519 0.105746 5.82 <1e-05 0.395887 0.834494
Subject: 9 0.674252 0.105746 6.38 <1e-05 0.454948 0.893555
Subject: 10 0.843855 0.105746 7.98 <1e-07 0.624551 1.06316
Subject: 11 0.32713 0.105746 3.09 0.0053 0.107826 0.546434
Subject: 12 0.556382 0.105746 5.26 <1e-04 0.337078 0.775686
Subject: 13 -0.277867 0.105746 -2.63 0.0154 -0.49717 -0.0585627
Subject: 14 -0.690512 0.105746 -6.53 <1e-05 -0.909816 -0.471208
Subject: 15 0.0 NaN NaN NaN NaN NaN
Subject: 16 -0.0551107 0.105746 -0.52 0.6075 -0.274415 0.164193
Subject: 17 -0.324189 0.105746 -3.07 0.0057 -0.543493 -0.104885
Subject: 18 -0.237008 0.105746 -2.24 0.0354 -0.456312 -0.017704

```

---

<div class="post-metadata">

**Author:** ![nalimilan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nalimilan/32/147_2.png) [@nalimilan](https://discourse.julialang.org/u/nalimilan)\
**Post date:** [July 25, 2022, 8:18am UTC](https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733/2 "2022-07-25T08:18:54Z")

</div>

I don’t really understand what you propose. Could you give more details? What kind of operator would you suggest?

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [July 25, 2022, 9:04am UTC](https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733/3 "2022-07-25T09:04:37Z")

</div>

For example, we have data and code bellow. Result will be simmilar to posted above. `Subject: 1` not included in model matrix (this is a base level), but `Subject: 15` included and it make collinearity. So, if make DummyCoding separetly for subjects in TR sequence and in RT sequence there will be no problems.

Code:

```julia
Using DataFrames, CSV, GLM
dfrds = CSV.File("rds1.csv", types = Dict(:Var => Float64, :Subject => String, :Period => String, :Sequence => String, :Formulation => String )) |> DataFrame
olscmax = lm(@formula(log(Var) ~ Sequence+Period+Formulation+Subject), dfrds; dropcollinear = true,)

```

Data:

```julia
Subject	Sequence	Period	Formulation	Var
3	TR	1	T	225.95
1	RT	1	R	181.09
2	RT	1	R	114.48
4	RT	1	R	176.91
5	TR	1	T	147.01
6	TR	1	T	97.53
7	RT	1	R	146.60
8	TR	1	T	45.58
9	RT	1	R	109.20
10	RT	1	R	125.61
11	TR	1	T	92.26
12	RT	1	R	237.95
13	TR	1	T	145.46
14	TR	1	T	179.96
15	TR	1	T	173.86
16	RT	1	R	144.00
17	RT	1	R	185.10
18	TR	1	T	117.99
1	RT	2	T	210.14
2	RT	2	T	98.72
3	TR	2	R	241.09
4	RT	2	T	186.65
5	TR	2	R	139.56
6	TR	2	R	124.77
7	RT	2	T	137.62
8	TR	2	R	57.71
9	RT	2	T	139.36
10	RT	2	T	120.43
11	TR	2	R	116.10
12	RT	2	T	228.63
13	TR	2	R	165.09
14	TR	2	R	181.09
15	TR	2	R	206.66
16	RT	2	T	143.25
17	RT	2	T	192.22
18	TR	2	R	125.50

```

---

<div class="post-metadata">

**Author:** ![Eric](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eric/32/25565_2.png) [@Eric](https://discourse.julialang.org/u/Eric)\
**Post date:** [July 25, 2022, 4:23pm UTC](https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733/4 "2022-07-25T16:23:54Z")

</div>

I tried in with [LinearRegressionKit.jl](https://github.com/ericqu/LinearRegressionKit.jl).

```julia
using LinearRegressionKit, StatsModels, DataFrames, CSV

dfrds = CSV.File("rds1.csv", types = Dict(:Var => Float64, :Subject => String, :Period => String, :Sequence => String, :Formulation => String )) |> DataFrame

olscmax = regress(@formula(log(Var) ~ Sequence+Period+Formulation+Subject), dfrds, req_stats=["default", "vif"])

```

and I get ok (as you will see, there are some warnings because some of your categories have no observation, which is probably because you have extracted only a part of the dataset). I also added the VIF to have something about the potential collinearity. I haven’t looked in much detail, but the lack of observations for some categories probably gives weird results in GLM. I would recommend using [FreqTables.jl](https://github.com/nalimilan/FreqTables.jl) from @nalimilan to get a more precision where the problem lies.

```julia
Model definition: :(log(Var)) ~ 1 + Sequence + Period + Formulation + Subject
Used observations: 23
Model statistics:
  R²: 0.991697 Adjusted R²: 0.908672
  MSE: 0.013198 RMSE: 0.114883
  σ̂²: 0.013198
Confidence interval: 95%

Coefficients statistics:
 Terms ╲ Stats │ Coefs Std err t Pr(>|t|) low ci high ci VIF
───────────────┼──────────────────────────────────────────────────────────────────────────────────────────
(Intercept) │ 5.26434 0.0938013 56.1223 0.000317339 4.86075 5.66793 0.0
Sequence: TR │ -8.2708 3.85482e6 -2.14557e-6 0.999998 -1.6586e7 1.6586e7 6.60528e15
Period: 2 │ 0.0122562 0.0741564 0.165275 0.883923 -0.306813 0.331326 1.63043
Formulation: T │ 0.00583091 0.0741564 0.0786299 0.944486 -0.313238 0.3249 2.3913
Subject: 10 │ -0.431159 0.148313 -2.90709 0.100759 -1.0693 0.20698 1.5942
Subject: 11 │ 7.52524 3.85482e6 1.95216e-6 0.999999 -1.6586e7 1.6586e7 1.10088e15
Subject: 12 │ 0.20772 0.148313 1.40055 0.296333 -0.430419 0.845859 1.5942
Subject: 13 │ 7.98053 3.85482e6 2.07027e-6 0.999999 -1.6586e7 1.6586e7 1.10088e15
Subject: 14 │ 8.19337 3.85482e6 2.12549e-6 0.999998 -1.6586e7 1.6586e7 1.10088e15
Subject: 15 │ 8.15888 3.85482e6 2.11654e-6 0.999999 -1.6586e7 1.6586e7 1.10088e15
Subject: 16 │ -0.294527 0.148313 -1.98585 0.185443 -0.932666 0.343612 1.5942
Subject: 17 │ -0.0434443 0.148313 -0.292923 0.797177 -0.681583 0.594694 1.5942
Subject: 18 │ 7.77123 3.85482e6 2.01598e-6 0.999999 -1.6586e7 1.6586e7 1.10088e15
Subject: 2 │ -0.60704 0.114883 -5.284 0.0339997 -1.10134 -0.11274 1.82609
Subject: 3 │ 8.45016 3.85482e6 2.1921e-6 0.999998 -1.6586e7 1.6586e7 2.10168e15
Subject: 4 │ -0.0709459 0.114883 -0.617551 0.599816 -0.565246 0.423354 1.82609
Subject: 5 │ 7.96192 3.85482e6 2.06544e-6 0.999999 -1.6586e7 1.6586e7 2.10168e15
Subject: 6 │ 7.58079 3.85482e6 1.96657e-6 0.999999 -1.6586e7 1.6586e7 1.10088e15
Subject: 7 │ -0.276633 0.148313 -1.8652 0.203152 -0.914771 0.361506 1.5942
Subject: 8 │ 6.8201 3.85482e6 1.76924e-6 0.999999 -1.6586e7 1.6586e7 1.10088e15
Subject: 9 │ -0.571159 0.148313 -3.85104 0.0612943 -1.2093 0.0669793 1.5942

```

Thought if you want to be specific about which category is the base you can specify it with the contrast (same format as in GLM). you can also find an example here: [regression example with contrasts](https://ericqu.github.io/LinearRegressionKit.jl/stable/multi_tutorial/)

```julia
olscmax = regress(@formula(log(Var) ~ Sequence+Period+Formulation+Subject), 
            dfrds, req_stats=["default", "vif"], 
            contrasts= Dict(:Subject => DummyCoding(base="1") ))

```

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [July 25, 2022, 7:28pm UTC](https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733/5 "2022-07-25T19:28:08Z")

</div>

Hi! Thanks a lot for [LinearRegressionKit.jl](https://github.com/ericqu/LinearRegressionKit.jl), it works perfect! 👌  
But this problem is not because data is incomplete. It is complete. This happens because we have in same time `Sequence` factor and `Subject`, and `Subject` is nested with `Sequence` - One level of `Subject` corresponds to only one level of `Sequence`. In some software (IBM SPSS for example) for this case used syntax like this: `/DESIGN=Sequence Period Formulation Subject(Sequence)` - this means that even `Subject` have same levels they coded separately for each `Sequence`, and also this means that DummyCoding (for n-1 factor levels) used for each `Sequence` (so if we have 2 sequences - we have n - 2 subject levels, if 3 sequences we have n - 3 subject levels ets.).

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [July 25, 2022, 10:08pm UTC](https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733/6 "2022-07-25T22:08:31Z")

</div>

So, one more example with data above:

```julia
Using DataFrames, CSV, GLM
model = lm(@formula(log(Var) ~ Sequence+Period+Formulation+Subject), dfrds; dropcollinear = true,)
using MixedAnova
glm_init()
MixedAnova.anova(model, type = 3)

```

ANOVA results:

```julia
Analysis of Variance

Type 3 test / F test

:(log(Var)) ~ 1 + Sequence + Period + Formulation + Subject

Table:
───────────────────────────────────────────────────────────
             DOF Exp.SS Mean Square F value Pr(>|F|)
───────────────────────────────────────────────────────────
(Intercept) 1 49.86 49.86 7795.7027 <1e-22
Sequence 1 1.7847 1.7847 279.0437 <1e-10
Period 1 0.0453 0.0453 7.0905 0.0170
Formulation 1 0.0228 0.0228 3.5725 0.0770
Subject 17
(Residuals) 16 0.1023 0.0064
───────────────────────────────────────────────────────────

```

We see: wrong DOF for `Subject`, absence SS and MS for `Subject`, wrong F value for `Sequence` (because Sequence is inter-subject factor and F should be = MSsequence / MSsubject), not residual. It is not problem of coding but it can help packages like MixedAnova.jl correctly handle this data structure.

And the question: do we need in StatsModels some operator like `var ~ factor1 + factor2 in factor1` to handle this models correctly?

---

<div class="post-metadata">

**Author:** ![Eric](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eric/32/25565_2.png) [@Eric](https://discourse.julialang.org/u/Eric)\
**Post date:** [July 25, 2022, 10:57pm UTC](https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733/7 "2022-07-25T22:57:32Z")

</div>

Thanks for providing more details.  
If I understand correctly, you’re looking for nested effects.  
I believe they are defined in `MixedModels.jl`; there is a [documentation here](https://juliastats.org/MixedModels.jl/stable/constructors/#Examples-of-linear-mixed-effects-model-fits) I believe they use the `/` as it builds on the syntax from R.

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [July 26, 2022, 12:04pm UTC](https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733/8 "2022-07-26T12:04:29Z")

</div>

Yep! I use it when possible MixedModels, Metida, GLM, LsqFit, FixedEffectsModels, ets… But sometime by statistical plan I should use ANOVA. And there is no problem for me - for some task I can use R or SPSS or write self-implemented anova. But even in self implemented tasks I found that I can’t make full-rank model matrix with StatsModels for some common cases. So maybe additional functionality for StatsModels help developers to make appropriate implementation of ANOVA…

---

<div class="post-metadata">

**Author:** ![palday](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palday/32/12640_2.png) [@palday](https://discourse.julialang.org/u/palday)\
**Post date:** [July 26, 2022, 3:27pm UTC](https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733/9 "2022-07-26T15:27:03Z")

</div>

@PharmCat have you seen [GitHub - kleinschmidt/RegressionFormulae.jl](https://github.com/kleinschmidt/RegressionFormulae.jl) ?

---

<div class="post-metadata">

**Author:** ![palday](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/palday/32/12640_2.png) [@palday](https://discourse.julialang.org/u/palday)\
**Post date:** [July 26, 2022, 3:32pm UTC](https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733/10 "2022-07-26T15:32:25Z")

</div>

I think part of the idea is to automatically drop columns corresponding to unestimable contrasts, e.g. when certain combinations of levels from different categorical predictors never co-occur – missing or empty cells in a quasi-factorial design. Currently, MixedModels.jl and GLM.jl both handle this case as part of more general pivoting / rank deficiency detection, but it would be good to support dropping clearly unestimable columns (e.g., constant and constant multiples of other columns) from the model matrix or even prevent the allocation of such columns. I’ve talked to @dave.f.kleinschmidt a bit about this and part of the problem at the moment IIRC is that this falls in awkwardly between `apply_schema` (schema-time) and `modelcols` (data-time) operations.

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [July 26, 2022, 4:02pm UTC](https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733/11 "2022-07-26T16:02:00Z")

</div>

Yes! Thank you. That’s exactly what i wanted to say. 😂 Second part of suggestion - decide - do we need operator in StatModels or in [RegressionFormulae](http://www.davekleinschmidt.com/RegressionFormulae.jl/dev/#RegressionFormulae.jl) that make model constructing more clear for nested factors. I look at RegressionFormulae it gives `^` and `/` - and it just expand term to others. Idea is make some operator for example `∈` with following mechanism:

model `var~ sequence + subject ∈ sequence`, contrasts = Dict(:subject =\> DummyCoding(base = Dict(“A” =\> 1, “B” =\> 4))))

data:

```julia
var sequence subject
. A 1
. A 2
. A 3
. B 4
. B 5
. B 6

```

should give model matrix:

```julia
1 0 0 0 0 0
1 0 1 0 0 0
1 0 0 1 0 0
1 1 0 0 0 0
1 1 0 0 1 0
1 1 0 0 0 1

```

that is full-ranked for coefs: int, sequence: B, subject(sequence: A): 2, subject(sequence: A): 3, subject(sequence: B): 5, subject(sequence: B): 6

---

<div class="post-metadata">

**Author:** ![dave.f.kleinschmidt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dave.f.kleinschmidt/32/55_2.png) [@dave.f.kleinschmidt](https://discourse.julialang.org/u/dave.f.kleinschmidt)\
**Post date:** [July 26, 2022, 4:59pm UTC](https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733/12 "2022-07-26T16:59:06Z")

</div>

To unpack this a bit more: currently, we have three stages between `@formula(...)` and actual matrices that get passed to GLM.jl etc.:

1. The formula expr is parsed into Terms (purely syntactic operations, nothing known about types of variables)
2. A schema is extracted from the data: categorical (i.e. non-numeric) variables are detected adn their unique values are extracted, one column at a time. This schema is then _applied_ to the parsed Terms, replacing `Term(:a)` with e.g. `CategoricalTerm(:a, ...)`.
3. These “concrete terms” can now be fed a table of data and will generate the numeric matrices

The issue is that, in order to support this kind of operator (create different contrasts given some kind of special syntax), we’d need to know about that syntax BEFORE the schema is computed; when the schema is applied, the full dataset is not available (by design). But special syntax is (at the moment) supported by operations that take place _during schema application_.

There are a number of reasons why this is not great, all with a similar shape to this: you want to somehow control the way you extract teh necessary invariants from the data to make the schema based on some special syntax/operators/model context. Some discussion here: [RFC: add an additional apply\_context stage · Issue #154 · JuliaStats/StatsModels.jl · GitHub](https://github.com/JuliaStats/StatsModels.jl/issues/154)

In the short term, I’m afraid the best solution is to create a new variable where you just `string` the pairs together; then the contrasts shoudl be computed correctly.

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [July 30, 2022, 11:07pm UTC](https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733/13 "2022-07-30T23:07:51Z")

</div>

Hi! Thank you for the great explanation!

I try to intercept `in` function term to make specific term with this function for better `nested ` contrast :

```julia
function StatsModels.apply_schema(
    t::FunctionTerm{typeof(∈)},
    sch::Schema,
    Mod::Type{<:RegressionModel},
    )
    length(t.args_parsed) == 2 ||
        throw(ArgumentError("malformed term: $t (Exactly two arguments required)"))

    #=
    How to get data here?
    =#
    #Make specific term NestedTerm(sym, contrasts, base)
end

```

But how to get access to data here? There are too many apply\_schema in StstsModels…

Or should it be done in `schema` stage? Or should I reload `concrete_term`?

StstsModels is a very confusing package - so I can’t to trace the way of `FunctionTerm{typeof(∈)}`

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [August 1, 2022, 9:25am UTC](https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733/14 "2022-08-01T09:25:07Z")

</div>

> [@dave.f.kleinschmidt](#):
>
> Some discussion here: [RFC: add an additional apply\_context stage · Issue #154 · JuliaStats/StatsModels.jl · GitHub](https://github.com/JuliaStats/StatsModels.jl/issues/154)

Hi! Could you kindly clarify - what is the status of this issue?

---

<div class="post-metadata">

**Author:** ![dave.f.kleinschmidt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dave.f.kleinschmidt/32/55_2.png) [@dave.f.kleinschmidt](https://discourse.julialang.org/u/dave.f.kleinschmidt)\
**Post date:** [August 1, 2022, 8:38pm UTC](https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733/15 "2022-08-01T20:38:17Z")

</div>

You’re lookin’ at it (in the issue 😉 ). It’s a pretty disruptive change so it needs some serious time to sit down and get it right and there’s a number of higher priority things. Plus I write julia all day long at my day job now so there’s not as much bandwidth to push this stuff forward as much as I’d like!

---

<div class="post-metadata">

**Author:** ![dave.f.kleinschmidt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dave.f.kleinschmidt/32/55_2.png) [@dave.f.kleinschmidt](https://discourse.julialang.org/u/dave.f.kleinschmidt)\
**Post date:** [August 1, 2022, 8:45pm UTC](https://discourse.julialang.org/t/should-model-matrix-for-nested-factors-be-full-rank/84733/16 "2022-08-01T20:45:00Z")

</div>

RE “how to get data here”: you can’t, and that’s by design. `apply_schema` does just that: applies teh information _in the schema_ to transform terms into different terms. If the information you need isn’t present in the schema, you can’t do what you’re trying to do I’m afraid.

The reason that constraint is important (IMO) is that we want to support streaming/online/out of core stuff too, or at least not _prevent_ support for it in the future. So lots of stuff that seems reasonable (like having access during `apply_schema` to the full dataset, which are usually pretty small data frames in teh grand scheme of things) are not allowed.

It’s an analogous situation to the splines situation ([Safe predictions using `@formula` and regression splines · Issue #168 · JuliaStats/StatsModels.jl · GitHub](https://github.com/JuliaStats/StatsModels.jl/issues/168)).

If we had a separate “apply context” stage, you could hook in there to create your `NestedTerm`, and then we could specialize schema extraction on these custom term types. I do think that’s a better design overall but like I said above, don’t really have the bandwidth to really think through how to make it work without adding TOO much extra complexity (which, as you already pointed out, already makes StatsModels.jl an intimidating package for many developers)