While several packages with functions that plot marginal effects already exist (such as interplot), these functions usually have limited applicability. For example, interplot does not allow the user to specify robust or clustered standard errors, nor does it plot interaction effects for censReg or robust regression objects produced using the robustbase package.

To solve this problem, I created a function called meplot(). This function not only allows the user to use custom variance-covariance matrices, it also accepts all sorts of model objects that existing options, like interplot, cannot accommodate.

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, estimate a regression model with an interaction term. .

mod <- lm(mpg ~ hp + wt + hp*wt, data=mtcars) # Using the mtcars dataset.

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

Then use the above to plot marginal effects.

library(ggplot2)
ggplot(data=NULL,aes(x=z0, y=dy.dx)) +
    labs(x="Weight",y="Marginal Effect of Horsepower",title="Example Marginal Effect Plot") +
    geom_line(aes(z0, dy.dx),size = 1) +
    geom_line(aes(z0, lwr), size = 1, 
              linetype = 2) +
    geom_line(aes(z0, upr), size = 1, 
              linetype = 2) +
    geom_hline(yintercept=0)

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 you want to plot a lot of marginal effects.

Enter meplot(). meplot() will do all the above math and plot marginal effects, all while requiring little effort on the part of the user. Even better, it allows the user to specify a modified variance-variance covariance matrix (which could be done “by hand” using the above code as a template).

First, the function itself:

# meplot()
meplot <- function(model,var1,var2,int,vcov,ci=.95,
                   xlab=var2,ylab=paste("Marginal Effect of",var1),
                   main="Marginal Effect Plot",
                   me_lty=1,me_lwd=1,me_col="black",
                   ci_lty=1,ci_lwd=.5,ci_col="black",
                   yint_lty=2,yint_lwd=1,yint_col="black"){
  require(ggplot2)
  alpha <- 1-ci
  z <- qnorm(1-alpha/2)
  beta.hat <- coef(model)
  cov <- vcov
  z0 <- seq(min(model.frame(model)[,var2],na.rm=T),max(model.frame(model)[,var2],na.rm=T),length.out=1000)
  dy.dx <- beta.hat[var1] + beta.hat[int]*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
  ggplot(data=NULL,aes(x=z0, y=dy.dx)) +
    labs(x=xlab,y=ylab,title=main) +
    geom_line(aes(z0, dy.dx),size = me_lwd, 
              linetype = me_lty, 
              color = me_col) +
    geom_line(aes(z0, lwr), size = ci_lwd, 
              linetype = ci_lty, 
              color = ci_col) +
    geom_line(aes(z0, upr), size = ci_lwd, 
              linetype = ci_lty, 
              color = ci_col) +
    geom_hline(yintercept=0,linetype=yint_lty,
               size=yint_lwd,
               color=yint_col)
}

The above function contains all the math and plotting parameters I used above, and requires only five user commands:

  1. A model object model (e.g. model=mod).
  2. The variable for which marginal effects are to be plotted var1 (e.g. var1="hp").
  3. The variable on which the marginal effects of the first is conditional var2 (e.g. var1="wt").
  4. The name of the interaction term, int (e.g. int="hp:wt").
  5. A variance-covariance object, vcov (e.g. vcov=robcov1).

By 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:

Below, I utilize the mod object I created above to demonstrate how meplot() works.

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

Also, because meplot() relies on ggplot2, any ggplot2 functions can be used to customize your marginal effect plot using the standard + geom_point(), for example.

meplot(model=mod,var1="hp",var2="wt",int="hp:wt",vcov=vcov(mod),
       xlab="Weight",ylab="Marginal Effect of Horesepower",
       main="A Better Marginal Effect Plot") + 
  theme_classic() +                        # Use the classic theme.
  theme(text=element_text(family="serif")) # Change the font.

Say that you need to use robust and/or clustered standard errors. meplot() can accommodate this. Simply specify a new variance-covariance matrix with something like the vcovHC() function in the sandwhich library.

library(sandwich)
robcov <- vcovHC(mod,type="HC3")  
meplot(model=mod,var1="hp",var2="wt",int="hp:wt",vcov=robcov)

That’s meplot()!