# How to find marginal effect of restricted cubic spline

I'm trying to figure out how to find the marginal effect of an interaction term from a restricted cubic spline in a non-linear model. The post Nonlinear effect in an interaction term is a good start on modeling the nonlinear effects and how to get plots, but does not address finding the marginal effect.

The package postrcspline in STATA has a function mfxrcspline which "displays the marginal effect of a restricted cubic spline," which is exactly what I am after. (See Figure 1 below)

R does not seem to offer this feature as conveniently ,so I'm trying to figure out how to get these same results.

As I understand it, suppose I have a multi-variable regression with restricted cubic splines and an interaction:

$$y = \beta_{0} + \beta_{1}x1 + \beta_{2} \mathcal{f}(x2) + \beta_{3} \mathcal{f}(x2) \cdot x1 + \epsilon$$

where $\mathcal{f}(x2)$ is a spline of the time-series (year)

The marginal effect of $\frac{\partial y}{\partial x1}$ is:

$$\frac{\partial y}{\partial x1} = \beta_{1} + \beta_{3} \mathcal{f}(x2)$$

where $\beta_{3}$ is the coefficient on the spline and $\mathcal{f}(x2)$ is a design matrix for each year in the regression that causes the slope to change for each $y$.

To say in words, I would like to find the marginal effect of $y$ for each year $x2$ in the spline given $\beta_{3}$.

In other words, it shows for each value of the spline variable how much the expected value of your explained variable changes for a unit change in the spline variable. It is the first derivative of the curve.

This appears to be simple matrix multiplication to plot the marginal effect, but I'm not sure how to statistically do this.

Here is a plot to illustrate what I'm after:

Figure 1: The left plot shows the results of the regression using a restricted cubic spline and the right provides the marginal effect--note the changes on the y-axis. Here is an R example to demonstrate the nonlinear effect from the regression (left plot in Figure 1):

library(rms)
set.seed(5)
# Fit a complex model and approximate it with a simple one
x1 <- runif(200)
x2 <- runif(200)
y <- x1 + x2 + rnorm(200)
f <- ols(y ~ x1 + rcs(x2,4)  + rcs(x2,4)*x1) 