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