I am trying to replicate some Stata code that uses average marginal effects to interpret interaction effects in R.

Suppose I am using the mtcars dataset to estimate the average miles/gallon mpg of a car as a function of gross horsepower hp that is allowed to vary by the number of cylinders cyl, like below:

mod1 <- lm(mpg ~ hp:as.factor(cyl), data=mtcars)

In Stata, I could calculate the average marginal effects of horsepower at different values of the discrete covariate cyl as follows:

margins, dydx(hp) at(cyl == 4) at(cyl == 6) at(cyl == 8)

I am wondering how I can do this in R, using the marginaleffects package.

I was thinking about using the avg_slopes function with the by statement (as show below), but I’m not sure whether the by statement matches the at statement of the Stata code. I don’t really understand the documentation of the by statement well.

library(marginaleffects)

avg_slopes(mod1, variables=c("hp"), by=c("cyl"))

**Output**
 Term    Contrast cyl Estimate Std. Error      z Pr(>|z|)    S   2.5 %   97.5 %
   hp mean(dY/dX)   4  0.00835     0.0347  0.241   0.8097  0.3 -0.0596  0.07633
   hp mean(dY/dX)   6 -0.04492     0.0248 -1.810   0.0704  3.8 -0.0936  0.00373
   hp mean(dY/dX)   8 -0.04736     0.0135 -3.496   <0.001 11.0 -0.0739 -0.02081