While several packages with functions that plot marginal effects already exist (such as interplot), these methods usually have limited applicability to only a certain set of regression objects. While plotting marginal effects can be done “by hand” in R rather easily, it would be best to avoid having to do this every time one runs into an issue using functions such as interplot() when one has to utilize models such as the lmrob() function in the robustbase package, which does not work with interplot(). To get around this issue, I detail below the meplot() function, which will take care of all the heavy lifting involved in creating a marginal effect plot, requiring the user to supply only the regression model object and the two variables that have been interacted in said model.

Doing It the Hard Way

Marginal effects (ME) are typically computed for multiple regression models where an interaction term is included, such as follows:

\[y_i=\beta_0+\beta_1x_i+\beta_2z_i + \beta_3(x_i \times z_i) + \epsilon_i\]

Once this model is estimated, the instantaneous effect of \(x_i\) on \(y_i\) as \(z_i\) varies (or vice versa) is calculated

\[\frac{\partial y}{\partial x} = \hat{\beta_x} + \hat{\beta}_{xz}z_i\]

along with standard errors

\[SE(\partial y/\partial x) = \sqrt{cov(\beta_x, \beta_x) + z_i^2cov(\beta_{xz}, \beta_{xz}) + 2z_icov(\beta_x, \beta_{xz})}\]

where the upper and lower bounds are estimated as follows:

\[\frac{\partial y}{\partial x} \pm t^*SE(\partial y/\partial x)\]

Doing the math in R is fairly straightforward. First, I estimate a regression model with an interaction term. Below, I utilize data from mtcars (one of the toy datasets in R), and I estimate a model where miles per gallon for a set of cars (mpg) is a function of horsepower (hp) and weight (wt).

mod <- lm(mpg ~ hp + wt + hp*wt, data=mtcars)

I then extract model estimates and the covariance matrix to compute marginal effects.

beta.hat <- coef(mod) 
cov <- vcov(mod)
z0 <- seq(min(mtcars$wt), max(mtcars$wt), length.out = 1000)
dy.dx <- beta.hat["hp"] + beta.hat["hp:wt"]*z0
se.dy.dx <- sqrt(cov["hp", "hp"] + z0^2*cov["hp:wt", "hp:wt"] + 2*z0*cov["hp", "hp:wt"])
upr <- dy.dx + 1.96*se.dy.dx
lwr <- dy.dx - 1.96*se.dy.dx

I then use the above to plot marginal effects.

par(family="serif",bty="l",mar=c(5,5.5,2,2))
plot(x=z0, y=dy.dx,type="n",xlim=c(min(z0),max(z0)),
     ylim=c(min(lwr),max(upr)),
     xlab = "Weight",
     ylab = expression(frac(partialdiff*paste("MPG"),
                            partialdiff*paste("Horsepower"))),
     main="Marginal Effect of Horsepower on MPG as Weight Varies")
lines(z0, dy.dx, lwd = 3)
lines(z0, lwr)
lines(z0, upr)
abline(h=0,lty=2)

Doing It the Easy Way

While doing all the above works just fine, “just fine” requires a lot of lines of code, which can get quite tedious if I want to plot a lot of marginal effects. Thankfully, I can do better than just fine.

Enter meplot(). meplot() will do all the above math and plot the marginal effects, all while requiring little effort on the part of the user. Let’s see it in action.

First, the function itself, with all its dirty innards:

# ------------------------------------
# Create Marginal Effect Plot Function
# ------------------------------------
meplot <- function(model,var1,var2,ci=.95,
                   xlab=var2,ylab=paste("Marginal Effect of",var1),
                   main="Marginal Effect Plot",
                   me_lty=1,me_lwd=3,me_col="black",
                   ci_lty=1,ci_lwd=1,ci_col="black",
                   yint_lty=2,yint_lwd=1,yint_col="black"){
  alpha <- 1-ci
  z <- qnorm(1-alpha/2)
  beta.hat <- coef(model)
  cov <- vcov(model)
  z0 <- seq(min(model$model[,var2],na.rm=T),max(model$model[,var2],na.rm=T),length.out=1000)
  dy.dx <- beta.hat[var1] + beta.hat[length(beta.hat)]*z0
  se.dy.dx <- sqrt(cov[var1,var1] + z0^2*cov[nrow(cov),ncol(cov)] + 2*z0*cov[var1,ncol(cov)])
  upr <- dy.dx + z*se.dy.dx
  lwr <- dy.dx - z*se.dy.dx
  plot(x=z0, y=dy.dx,type="n",xlim=c(min(z0),max(z0)),
       ylim=c(min(lwr),max(upr)),
       xlab = xlab,
       ylab = ylab,
       main = main)
  lines(z0, dy.dx, lwd = me_lwd, lty = me_lty, col = me_col)
  lines(z0, lwr, lwd = ci_lwd, lty = ci_lty, col = ci_col)
  lines(z0, upr, lwd = ci_lwd, lty = ci_lty, col = ci_col)
  abline(h=0,lty=yint_lty,lwd=yint_lwd,col=yint_col)
}

The above function essentially contains all the math and plotting parameters I used above, and requires only three user commands: a model object model, the variable for which marginal effects are to be plotted var1, and the variable on which the marginal effects of the first is conditional var2. Be default, meplot() computes 95% confidence intervals; however the user can change this by specifying a different value for ci. X and Y axis labels, along with several other plotting parameters, can also be changed, including:

Let’s see it in action. Below, I utilize the mod object I created above.

meplot(model=mod,var1="hp",var2="wt")

Also, meplot() will work with any base R plotting parameters specifiable with the par() function, which is good to know if the default parameters aren’t quite your style.

par(family="serif",bty="l",mar=c(5,5.5,2,2))
meplot(model=mod,var1="hp",var2="wt",xlab = "Weight", ci=.8,
     ylab = expression(frac(partialdiff*paste("MPG"),
                            partialdiff*paste("Horsepower"))),
     main="Marginal Effect of Horsepower on MPG as Weight Varies",
     ci_lty=3,yint_lty=1,yint_col="grey")

That’s meplot() in a nutshell!