library(marginaleffects)
dat <- transform(mtcars, y = mpg, x = wt)
mod <- lm(y ~ x, data = dat)
pct_change_x <- mean(dat$x) / 100
p0 <- predict(mod)
p1 <- predict(mod, newdata = transform(dat, x = x + pct_change_x))
mean(p1 - p0)
#> [1] -0.171945
s <- avg_slopes(mod, slope = "dyex")
s$estimate
#> [1] -17.194516 Elasticity
Slopes and elasticities can only be calculated for continuous numeric variables. The slopes() functions will automatically revert to comparisons() for binary or categorical variables.
In some contexts, it is useful to interpret the results of a regression model in terms of elasticity or semi-elasticity. One strategy to achieve that is to estimate a log-log or a semilog model, where the left and/or right-hand side variables are logged. Another approach is to note that \(\frac{\partial ln(x)}{\partial x}=\frac{1}{x}\), and to post-process the marginal effects to transform them into elasticities or semi-elasticities.
For example, say we estimate a linear model of this form:
\[y = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \varepsilon\]
Let \(\hat{y}\) be the adjusted prediction made by the model for some combination of covariates \(x_1\) and \(x_2\). The slope with respect to \(x_1\) (or “marginal effect”) is:
\[\frac{\partial \hat{y}}{\partial x_1}\]
We can estimate the slope (“dydx”), semi-elasticities (“dyex” and “eydx”), or elasticity (“eyex”) with respect to \(x_1\) as follows:
\[\eta_1=\frac{\partial \hat{y}}{\partial x_1}\] \[\eta_2=\frac{\partial \hat{y}}{\partial x_1}\cdot x_1\] \[\eta_3=\frac{\partial \hat{y}}{\partial x_1}\cdot \frac{1}{\hat{y}}\] \[\eta_4=\frac{\partial \hat{y}}{\partial x_1}\cdot \frac{x_1}{\hat{y}}\]
with interpretations roughly as follows:
- “dydx”: Increasing the \(x_1\) variable by 1 unit is associated with a change of \(\eta_1\) units in \(y\).
- “dyex”: Increasing the \(x_1\) variable by 100% of its baseline value is associated with a change of \(\eta_2\) units in \(y\).
- “eydx”: Increasing the \(x_1\) variable by 1 unit is associated with a change of \(\eta_3\)% of \(y\)’s baseline value.
- “eyex”: Increasing \(x_1\) by 1% of its baseline value is associated with a change of \(\eta_4\)% in the baseline value of \(y\).
With the marginaleffects package, these quantities are easy to compute with the slopes() function and the slope argument.
import numpy as np
import statsmodels.formula.api as smf
from marginaleffects import avg_slopes
from statsmodels.datasets import get_rdataset
mtcars = get_rdataset("mtcars").data
dat = mtcars.assign(
y=mtcars["mpg"],
x=mtcars["wt"]
)
mod = smf.ols("y ~ x", data=dat).fit()
pct_change_x = dat["x"].mean() / 100
p0 = mod.predict(dat)
p1 = mod.predict(dat.assign(x=dat["x"] + pct_change_x))
print(np.mean(p1 - p0))
#> -0.17194501167342113
s_dyex = avg_slopes(mod, slope="dyex")
print(float(s_dyex["estimate"][0]))
#> -17.194501167309795Increasing the x variable by 100% of its baseline value is associated with a change of \(\eta_1=-17.1945012\) units in \(y\). Notice that the slopes() estimates are scaled differently than in the manual computation.
yhat = mod.predict(dat)
p0 = mod.predict(dat)
p1 = mod.predict(dat.assign(x=dat["x"] + 1))
print(np.mean((p1 - p0) / yhat))
#> -0.2921699751468847
s_eydx = avg_slopes(mod, slope="eydx")
print(float(s_eydx["estimate"][0]))
#> -0.2921699751463526Increasing the x variable by 1 unit is associated with a change of \(\eta_2=-0.29217\%\) of \(y\)’s baseline value.
pct_change_x <- 1 / dat$x
pct_change_y <- (p1 - p0) / yhat
mean(pct_change_y / pct_change_x)
#> [1] -1.038292
s <- avg_slopes(mod, slope = "eyex")
s$estimate
#> [1] -1.038292pct_change_x = 1 / dat["x"]
pct_change_y = (p1 - p0) / yhat
print(np.mean(pct_change_y / pct_change_x))
#> -1.0382921374790237
s_eyex = avg_slopes(mod, slope="eyex")
print(float(s_eyex["estimate"][0]))
#> -1.0382921374770178Increasing x by 1% of its baseline value is associated with a change of \(\eta_3=-1.0382921\%\) in the baseline value of \(y\).