How to use delta method for standard errors of marginal effects?

Question

I am interested in better understanding the delta method for approximating the standard errors of the average marginal effects of a regression model that includes an interaction term. I've looked at related questions under delta-method but none have provided quite what I'm looking for.

Consider the following example data as a motivating example:

set.seed(1)
x1 <- rnorm(100)
x2 <- rbinom(100,1,.5)
y <- x1 + x2 + x1*x2 + rnorm(100)
m <- lm(y ~ x1*x2)

I am interested in the average marginal effects (AMEs) of x1 and x2. To calculate these, I simply do the following:

cf <- summary(m)$coef
me_x1 <- cf['x1',1] + cf['x1:x2',1]*x2 # MEs of x1 given x2
me_x2 <- cf['x2',1] + cf['x1:x2',1]*x1 # MEs of x2 given x1
mean(me_x1) # AME of x1
mean(me_x2) # AME of x2

But how do I use the delta method to calculate the standard errors of these AMEs?

I can calculate the SE for this particular interaction by hand:

v <- vcov(m)
sqrt(v['x1','x1'] + (mean(x2)^2)*v['x1:x2','x1:x2'] + 2*mean(x2)*v['x1','x1:x2'])

But I don't understand how to use the delta method.

Ideally, I'm looking for some guidance on how to think about (and code) the delta method for AMEs of any arbitrary regression model. For example, this question provides a formula for the SE for a particular interaction effect and this document from Matt Golder provide formulae for a variety of interactive models, but I want to better understand the general procedure for calculating SEs of AMEs rather than the formula for the SE of any particular AME.

Answer 1

The delta method simply says that if you can represent an auxiliary variable you can represent as a function of normally distributed random variables, that auxiliary variable is approximately normally distributed with variance corresponding to how much the auxiliary varies with respect to the normal variables. (EDIT: as pointed out by Alecos Papadopoulos the delta method can be stated more generally such that it does not require asymptotic normality.)

The easiest way to think of this is as a Taylor expansion, where the first term of a function is the mean, and the variance comes from the second order terms. Specifically, if $g$ is a function of parameter $\beta$ and $b$ is a consistent, normally distributed estimator for that parameter: $$ g(b) \approx g(\beta) + \nabla g(\beta)^\prime (b - \beta) $$ Since $\beta$ is a constant, and $b$ is a consistent estimator for $\beta$, we can then say: $$ \sqrt{n}\left(g(b)-g(\beta)\right)\,\xrightarrow{D}\,N\left(0, \nabla g(\beta)^\prime \cdot \Sigma_b \cdot \nabla g(\beta)\right) $$ In this case, $b$ is your OLS estimate, and $g$ is the AME. You can write this specific AME as: $$ g(b_1,b_2)=b_1+b_2 \text{ mean}(x_2) $$ if you took the gradient of this function (remember, a function of the coefficients not of $x_2$), it would be: $$ [1,\,\, \text{mean}(x_2)]^\prime $$ and the variance-covariance matrix for $b$ might be: $$ \left[ \begin{matrix} s_{11} & s_{12} \\ s_{12} & s_{22} \end{matrix}\right] $$ Plugging this into the variance formula and doing some matrix algebra gives you the same expression you wanted.

In general if you want to do this, you can explicitly code whatever $g$ you want into R as a function of all your coefficients and then use numDeriv to take the numerical gradient (otherwise you'd have to use computer algebra) of the function with respect to your parameters, at the parameters you estimated. Then you simply take the variance-covariance matrix and this numerical gradient and plug it into the formula and voila! Delta method.

ADDENDUM: In this specific case the R code would be:

v <- vcov(m)

# Define function of coefficients. Note all coefficients are included so it 
# will match dimensions of regression coefficients, this could be done more 
# elegantly in principle
g <- function(b){
    return(b[2] + b[4] * mean(x2))
}

require(numDeriv) # Load numerical derivative package

grad_g <-  jacobian(g, m$coef) # Jacobian gives dimensions, otherwise same as
                               # gradient 

sqrt(grad_g%*% v %*% t(grad_g)) # Should be exactly the same 

Note that it will always be preferable to get the exact gradient instead of the numerical gradient for this problem, as the exact gradient will have less computational error. The fact that $g$ is linear eliminates this problem, and for more complicated functions the exact gradient may not always be available.