# R utility for pulling out parameters for individual variable in lm when there are effect modifiers (interactions with categorical variables)

NOTE: I think the emmeans package may be what I'm looking for. Still welcome any input!

Suppose I have a regression model $$y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_1 x_2$$, where $$x_1$$ is continuous and $$x_2$$ has two levels: "A" and "B". MRE below.

From the lm summary, I can find the slope of y vs x for each factor level (call them groups A and B).

Group A: slope = $$\beta_1$$ = -0.034

Group B: slope = $$\beta_1 + \beta_3$$ = -0.034 + 0.05 = 0.016

I have some complex regression models with multiple continuous variables and effect modifiers and I would like to easily make a table of the slopes for the different subgroups. Rather than doing it by hand as above, is there a utility that extracts this information automatically?

I don't mind writing something myself, but I'd rather not reinvent the wheel when there might be a better one out there already.

library(tidyverse)
library(ggplot2)

set.seed(100)

b0 = 3
b1 = -0.034
b2 = -2
b3 = 0.05

x1 = 1:100
x2 = sample(c(0,1), 100, replace=T)

y = b0 + b1*x1 + b2*x2 + b3*x1*x2 + rnorm(100)

df = data.frame(Y=y, X1=x1, X2=ifelse(x2==0, 'A', 'B') %>% factor(levels=c('A', 'B')))

df %>%
ggplot(aes(x=X1, y=Y, color=X2)) +
geom_point() +
geom_smooth(method='lm')


df %>%
lm(y ~ x1*x2, .) %>%
summary()
#>
#> Call:
#> lm(formula = y ~ x1 * x2, data = .)
#>
#> Residuals:
#>      Min       1Q   Median       3Q      Max
#> -2.12027 -0.58294 -0.03404  0.58940  2.51982
#>
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)
#> (Intercept)  2.843017   0.319058   8.911 3.25e-14 ***
#> x1          -0.035285   0.005352  -6.592 2.34e-09 ***
#> x2          -1.775823   0.415229  -4.277 4.48e-05 ***
#> x1:x2        0.050924   0.007115   7.157 1.65e-10 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Residual standard error: 1.013 on 96 degrees of freedom
#> Multiple R-squared:  0.4292, Adjusted R-squared:  0.4114
#> F-statistic: 24.06 on 3 and 96 DF,  p-value: 1.071e-11


Created on 2020-03-19 by the reprex package (v0.3.0)

Yes, emmeans can do this.

> lm(y ~ x1 * x2) -> mod

> library(emmeans)
> emt <- emtrends(mod, "x2", var = "x1")
> emt
x2 x1.trend      SE df lower.CL upper.CL
0  -0.0353 0.00535 96 -0.04591  -0.0247
1   0.0156 0.00469 96  0.00633   0.0249

Confidence level used: 0.95

> pairs(emt)
contrast estimate      SE df t.ratio p.value
0 - 1     -0.0509 0.00711 96 -7.157  <.0001


That done, may I suggest using meaningful variable names? Who cares about x1 and x2? This is not a math class, it is a statistics forum. The meaning of what you are doing would be much more adequately explained using variable names like response, treatment, and dose.

• Thanks! Funny you should mention variable names. I am constantly drilling colleagues about self-commenting code. My variable names here were just to keep things super general for an example. Fear not--millitant missionary for legible code here! Mar 21, 2020 at 4:11
• Glad to hear it. Mar 21, 2020 at 13:19