2
$\begingroup$

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)

$\endgroup$

1 Answer 1

1
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ 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! $\endgroup$
    – abalter
    Mar 21, 2020 at 4:11
  • $\begingroup$ Glad to hear it. $\endgroup$
    – Russ Lenth
    Mar 21, 2020 at 13:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.