# Getting RegressionFormulae ^ (interaction) operator to work

**URL:** <https://discourse.julialang.org/t/getting-regressionformulae-interaction-operator-to-work/79216>\
**Category:** Statistics\
**Created:** [April 8, 2022, 9:10am UTC](https://discourse.julialang.org/t/getting-regressionformulae-interaction-operator-to-work/79216 "2022-04-08T09:10:35Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ross\_Boylan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ross_boylan/32/9210_2.png) [@Ross\_Boylan](https://discourse.julialang.org/u/Ross_Boylan)\
**Post date:** [April 8, 2022, 9:10am UTC](https://discourse.julialang.org/t/getting-regressionformulae-interaction-operator-to-work/79216/1 "2022-04-08T09:10:36Z")

</div>

My effort to use the `^` operator from [RegressionFormulae](https://github.com/kleinschmidt/RegressionFormulae.jl) doesn’t seem to have worked. What am I missing?

The docs just referenced say

> `(a+b+c)^2` generates `a + b + c + a&b + a&c + b&c` , but not `a&b&c`

Test code:

```julia
using CategoricalArrays
using StatsModels
using DataFrames
using RegressionFormulae

function follow(f::FormulaTerm, df::DataFrame)
    print("formula: ", f, "\n with terms ", terms(f))
    s = schema(f, df)
    print("\nschema: ", s)
    ts = apply_schema(f, s)
    print("apply_schema: ", ts, " with coefficient names", coefnames(ts))

    x = modelcols(ts.rhs, df)
    print("\nmodelcols: the the result has size ", size(x), " and type ", typeof(x))
    ## it's Matrix{Float64} and I see no labels on it. Nor does the type allow any.
    
end

df = DataFrame(a=[1, 2, 3], b=[0, 0, 1], c=[0.5, 2.5, 20], d=categorical(["xa", "xa", "xb"]))
follow(@formula(a~(a+b+c)^2), df)

```

Result:

```julia
formula: a ~ :((a + b + c) ^ 2)
  with terms Term[a, b, c]
schema: StatsModels.Schema with 3 entries:
  c => c
  b => b
  a => a
apply_schema: a ~ :((a + b + c) ^ 2) with coefficient names("a", "(a + b + c) ^ 2")
modelcols: the the result has size (3, 1) and type Matrix{Float64}

```

Expected result:  
After `apply_schema` coefficient names would reflect the expanded list of terms, and `modelcols` would have more than 1 column (like 7).

---

<div class="post-metadata">

**Author:** ![pdeffebach](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pdeffebach/32/10320_2.png) [@pdeffebach](https://discourse.julialang.org/u/pdeffebach)\
**Post date:** [April 8, 2022, 8:06pm UTC](https://discourse.julialang.org/t/getting-regressionformulae-interaction-operator-to-work/79216/2 "2022-04-08T20:06:16Z")

</div>

Interesting. A cursory glance at the package indicates this may be a bug. Please file an issue at RegressionFormulae.jl

---

<div class="post-metadata">

**Author:** ![Ross\_Boylan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ross_boylan/32/9210_2.png) [@Ross\_Boylan](https://discourse.julialang.org/u/Ross_Boylan)\
**Post date:** [April 9, 2022, 5:23am UTC](https://discourse.julialang.org/t/getting-regressionformulae-interaction-operator-to-work/79216/3 "2022-04-09T05:23:58Z")

</div>

Short version: I think `RegressionFormulae` requires contextual cues for the interaction operator to work. I didn’t provide them.

The test suite for `RegressionFormulae` for “powers of terms” has some crucial differences from what I did. It also has some broken tests covering more exotic cases than mine. I think the crucial thing is that they all get results by performing a `fit` on a special [dummy model](https://github.com/kleinschmidt/RegressionFormulae.jl/blob/main/test/dummymod.jl) defined for testing, and that model **is a subtype of `RegressionModel`**. The tests call `fit` on the model, which I think ends up invoking [this line](https://github.com/JuliaStats/StatsModels.jl/blob/bbd0b867acece8fb1051771ab6e20384def4f42b/src/statsmodel.jl#L78), which has a specialization for `RegressionModel`.

The code of that function invokes `ModelFrame` with the argument `model=RegressionModel`. I think as currently written `RegressionFormulae`s interpretation of `^` only operates when this [contextual cue](https://github.com/kleinschmidt/RegressionFormulae.jl/blob/934aae9c7cfe1e613bd8d091506a99ee3b0344bc/src/power.jl#L24) is available.

In contrast, my code invoked `schema`, `apply_schema`, and `modelcols` consistent with current recommendations for [`StatsModels`](https://juliastats.org/StatsModels.jl/latest/internals/); those methods are even in commented out code in the previous `fit` function. At least as I did it, they got no clues they were operating on a `RegressionModel`.

P.S. The one model column my test code produces is the result interpreting the formula very naively; it is the result of adding the 3 columns `a`, `b` and `c` and then squaring the result, element by element. So it is not only `^` but `+` that get interpreted as ordinary arithmetic, rather than special terms in the formula.

---

<div class="post-metadata">

**Author:** ![Ross\_Boylan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ross_boylan/32/9210_2.png) [@Ross\_Boylan](https://discourse.julialang.org/u/Ross_Boylan)\
**Post date:** [April 29, 2022, 6:41pm UTC](https://discourse.julialang.org/t/getting-regressionformulae-interaction-operator-to-work/79216/4 "2022-04-29T18:41:20Z")

</div>

If I change the invocation from

```julia
ts = apply_schema(f, s)

```

to

```julia
ts = apply_schema(f, s, RegressionModel)

```

(the last argument is the type that the modified handling of `^` applies to), things work as expected. The output is now

```julia
formula: a ~ :((a + b + c) ^ 2)
  with terms Term[a, b, c]
schema: StatsModels.Schema with 3 entries:
  b => b
  a => a
  c => c
apply_schema: a ~ 1 + a + b + c + a & b + a & c + b & c with coefficient names("a", ["(Intercept)", "a", "b", "c", "a & b", "a & c", "b & c"])
modelcols: the the result has size (3, 7) and type Matrix{Float64}
Model Matrix:
[1.0 1.0 0.0 0.5 0.0 0.5 0.0; 1.0 2.0 0.0 2.5 0.0 5.0 0.0; 1.0 3.0 1.0 20.0 3.0 60.0 20.0]

```

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [April 29, 2022, 7:20pm UTC](https://discourse.julialang.org/t/getting-regressionformulae-interaction-operator-to-work/79216/5 "2022-04-29T19:20:20Z")

</div>

Thanks for bringing this up. I raised an issue in this repo 2 years ago.  
If we want ppl to use Julia for serious stats work I think it is crucial to make these types of operations composable for stats-work throughout the entire Julia ecosystem:

> <https://github.com/kleinschmidt/RegressionFormulae.jl/issues/1>
>
> Hi @kleinschmidt,
> I think the Julia ecosystem would benefit from something like… this! 
> If we wanna do serious stats it should be easy to automatically generate all interactions (up order n) etc. 
> 
> Some things I find particularly useful in my other stats packages outside Julia: 
> 1. "i.x1" makes x1 into a \[factor variable\](https://www.stata.com/features/overview/factor-variables/) in a formula
> Suppose x1 takes the values: 1.2, 5, 6.4
> \`reg y x1\`: treats x1 as continuous & returns 1 coef (assuming no intercept)
> \`reg y i.x1\` creates 3 dummies for each level of x1 & returns 3 coefficients
> (if there is an intercept it randomly drops one level unless the user chooses which level to drop)
> 
> 2. \`i.x1#(c.x2 i.x3)\`
> Interacts all dummies of x1 w/ x2 (continuous)
> Interacts all dummies of x1 w/ all dummies of x3 
> 
> 3. Leads & Lags. Suppose D is at the state-year level. 
> \`L.D\`: creates a 1 year lag of D
> \`L(4).D\`: creates a 4 year lag of D
> \`F(4).D\`: creates a 4 year lead of D. $D\_{t+4}$
> \`reg y F(-1 0 1 2).D\`
> estimates: y\_t =b\_{-1} x\_{t-1} + b\_{0} x\_{t} +b\_{1} x\_{t+1} +b\_{2} x\_{t+2} 
> 
> If Julia is to be "\[as easy for statistics as R\](https://julialang.org/blog/2012/02/why-we-created-julia/)" these features should be in StatsModels.
> I'd love to help if I can.

@bkamins are there convenient-composable ways to do these things (raise a formula to a power, create interactions w/ factor variables, make arbitrary lead/lags) in the Julia ecosystem?  
-By composable I mean I can do the same thing whether I’m working w/ FixedEffects.jl or a time-series package…

---

<div class="post-metadata">

**Author:** ![bkamins](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bkamins/32/208538_2.png) [@bkamins](https://discourse.julialang.org/u/bkamins)\
**Post date:** [April 30, 2022, 6:35am UTC](https://discourse.julialang.org/t/getting-regressionformulae-interaction-operator-to-work/79216/6 "2022-04-30T06:35:10Z")

</div>

> [@Albert\_Zevelev](#):
>
> are there convenient-composable ways to do these things

I think we still need to work on this. Probably @kleinschmidt is the best person to ask about the current status.

---

<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:** [April 30, 2022, 8:53pm UTC](https://discourse.julialang.org/t/getting-regressionformulae-interaction-operator-to-work/79216/7 "2022-04-30T20:53:32Z")

</div>

`lead` and `lag` are supported by StatsModels: [Temporal variables and Time Series Terms · StatsModels.jl](https://juliastats.org/StatsModels.jl/stable/temporal_terms/). Interactions with factor (categorical) variables can be obtained using `a*b` or `a&b` just like with numeric variables. Do you have particular requirements that these options don’t meet?

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [May 1, 2022, 11:17pm UTC](https://discourse.julialang.org/t/getting-regressionformulae-interaction-operator-to-work/79216/8 "2022-05-01T23:17:55Z")

</div>

I went into some detail in the [issue](https://github.com/kleinschmidt/RegressionFormulae.jl/issues/1) I raised in the link above.  
There are certain basic operations that economists use daily in STATA/R that are not yet convenient in the Julia ecosystem:

One example:  
in STATA `reg y F(-1 0 1 2).D`  
estimates y\_{t} = \beta\_{-1} x\_{t-1} + \beta\_{0} x\_{t} + \beta\_{1} x\_{t+1} + \beta\_{2} x\_{t+2} +\varepsilon\_{t}  
I don’t know how to automatically do this (& the other things raised in the issue) in Julia

---

<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:** [May 13, 2022, 3:10pm UTC](https://discourse.julialang.org/t/getting-regressionformulae-interaction-operator-to-work/79216/9 "2022-05-13T15:10:06Z")

</div>

yeah, you got it 🙂 I think this is a docs issue if you had to go into the tests to figure it out. Specializing on the ‘context clues’ (model type we’re generating a formula for) is an important part of the mechanism that allows potentially competing claims to the same syntax.

we’d also consider changing RegressionFormulae.jl to use `Mod::Type{<:Any}` until someone complains.

also note that there are other situations where the behavior of `apply_schema` changes depending on the context, most notable the logic that promotes categorical terms to “full dummy” coding when necessary (only happens for `Mod::Type{<:StatisticalModel}` since it may or may not make sense to do that promotion in more ML-flavored situations).