# Running counterfactuals on a regression with interactions

**URL:** <https://discourse.julialang.org/t/running-counterfactuals-on-a-regression-with-interactions/56711>\
**Category:** New to Julia\
**Tags:** dataframes, glm\
**Created:** [March 8, 2021, 2:43am UTC](https://discourse.julialang.org/t/running-counterfactuals-on-a-regression-with-interactions/56711 "2021-03-08T02:43:55Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![samerb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samerb/32/9242_2.png) [@samerb](https://discourse.julialang.org/u/samerb)\
**Post date:** [March 8, 2021, 2:43am UTC](https://discourse.julialang.org/t/running-counterfactuals-on-a-regression-with-interactions/56711/1 "2021-03-08T02:43:55Z")

</div>

Let’s say I have a data frame with columns c1,c2,…,cn and I run a regression of c1 on c2, c3, c4 and c2.\*c4:

```julia
sm = GLM.lm( @formula(c1 ~ c2 + c3 +c4 +c2&c4), df)

```

If I were in Python, I could then do

```julia
formula = 'c1 ~ c2 + c3 +c4 +c2&c4'
import patsy
y, X = patsy.dmatrices(formula, df)

```

This would give me `X`, the design matrix of the regression–i.e. columns corresponding to c2, c3, c4, and c2.\*c4.  
How can I get this matrix in Julia? I could probably do this manually, but if I had a large number of regressors and interactions, it might get tricky, so I’d like a systematic approach.

The reason I’m asking is that I eventually want to do counterfactuals for what would happen if I change c2 to 0 or to 1. You could throw this in as a bonus, but I think once I get the first part, I can figure that out.

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [March 8, 2021, 2:53am UTC](https://discourse.julialang.org/t/running-counterfactuals-on-a-regression-with-interactions/56711/2 "2021-03-08T02:53:12Z")

</div>

I think you are looking for `modelmatrix(sm)`

---

<div class="post-metadata">

**Author:** ![samerb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samerb/32/9242_2.png) [@samerb](https://discourse.julialang.org/u/samerb)\
**Post date:** [March 8, 2021, 2:57am UTC](https://discourse.julialang.org/t/running-counterfactuals-on-a-regression-with-interactions/56711/4 "2021-03-08T02:57:58Z")

</div>

That doesn’t appear to be in the documentation. Just curious how you found it (not that important though)

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [March 8, 2021, 2:58am UTC](https://discourse.julialang.org/t/running-counterfactuals-on-a-regression-with-interactions/56711/5 "2021-03-08T02:58:58Z")

</div>

My brain 🙂

No idea when I found it.

It is in the StatsModels.jl docs. If something isn’t in GLM.jl docs you should look there.

[https://juliastats.org/StatsModels.jl/stable/api/#StatsBase.modelmatrix](https://juliastats.org/StatsModels.jl/stable/api/#StatsBase.modelmatrix)

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [March 8, 2021, 7:38am UTC](https://discourse.julialang.org/t/running-counterfactuals-on-a-regression-with-interactions/56711/6 "2021-03-08T07:38:21Z")

</div>

Why do you need the matrix though if you want to do counterfactuals? `predict` does this for you automatically:

```julia
julia> using DataFrames, GLM

julia> df = DataFrame(a = 1:10, b = 6:15);

julia> df.y = 0.5 * df.a .+ 0.8*df.b .+ df.a .* df.b;

julia> sm = lm(@formula(y ~ a*b), df)
StatsModels.TableRegressionModel{LinearModel{GLM.LmResp{Vector{Float64}}, GLM.DensePredChol{Float64, LinearAlgebra.CholeskyPivoted{Float64, Matrix{Float64}}}}, Matrix{Float64}}

y ~ 1 + a + b + a & b

Coefficients:
────────────────────────────────────────────────────────────────────────────────────
             Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95%
────────────────────────────────────────────────────────────────────────────────────
(Intercept) 0.0 NaN NaN NaN NaN NaN
a 0.5 4.8306e-13 1035067731159.83 <1e-81 0.5 0.5
b 0.8 1.23567e-13 6474215637663.89 <1e-87 0.8 0.8
a & b 1.0 2.28609e-14 43742897935693.22 <1e-93 1.0 1.0
────────────────────────────────────────────────────────────────────────────────────

julia> round.(predict(sm, DataFrame(a = 1:10, b = 0)))
10-element Vector{Float64}:
 0.0
 1.0
 1.0
 2.0
 2.0
 3.0
 3.0
 4.0
 4.0
 5.0

```