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 but none have provided quite what I'm looking for.

Consider the following example data as a motivating example:

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.


1 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.

  • 1
    $\begingroup$ Thanks for this very detailed answer. I think what was especially tripping me up was gradients with respect to the coefficients rather than the original variables. I really appreciate your help! $\endgroup$
    – Thomas
    Nov 4, 2014 at 6:15
  • 1
    $\begingroup$ It's equivalent for linear AMEs, when you take the average over the observations you just end up with the marginal effect at the mean. Otherwise you would really have to define g as the average of the marginal effects for each individual, and probably use the numerical gradient, I'm not sure that taking the SE for each would be quite the same. $\endgroup$
    – jayk
    Nov 4, 2014 at 11:56
  • 1
    $\begingroup$ That is AME and ME at the mean are equivalent for linear MEs. The SE won't I think be equivalent because the form for variance is quadratic, so the mean won't just pop out. I don't have good intuition for why the SE can't just be added up over observations, but I'm pretty sure it's true. $\endgroup$
    – jayk
    Nov 4, 2014 at 12:04
  • 2
    $\begingroup$ Note that the Delta Theorem does not require normality. $\endgroup$ Nov 5, 2014 at 2:26
  • 1
    $\begingroup$ @BerndWeiss I'm late to this, but this whole Q&A was a setup for my implementation thereof: cran.r-project.org/package=margins $\endgroup$
    – Thomas
    Aug 1, 2018 at 2:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.