¿Cómo utilizar el método delta para errores estándar de efectos marginales?

Estoy interesado en comprender mejor el método delta para aproximar los errores estándar de los efectos marginales promedio de un modelo de regresión que incluye un término de interacción. He analizado preguntas relacionadas con el método delta, pero ninguna me ha proporcionado lo que estoy buscando.

Considere los siguientes datos de ejemplo como un ejemplo motivador:

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

Estoy interesado en los efectos marginales promedio (AME) de x1y x2. Para calcular esto, simplemente hago lo siguiente:

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

Pero, ¿cómo uso el método delta para calcular los errores estándar de estos AME?

Puedo calcular el SE para esta interacción particular a mano:

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

Pero no entiendo cómo usar el método delta.

Idealmente, estoy buscando alguna guía sobre cómo pensar (y codificar) el método delta para AMEs de cualquier modelo de regresión arbitraria. Por ejemplo, esta pregunta proporciona una fórmula para el SE para un efecto de interacción particular y este documento de Matt Golder proporciona fórmulas para una variedad de modelos interactivos, pero quiero comprender mejor el procedimiento general para calcular los SE de AMEs en lugar de la fórmula para el SE de cualquier AME en particular.

El método delta simplemente dice que si puede representar una variable auxiliar, puede representarla en función de variables aleatorias normalmente distribuidas, esa variable auxiliar se distribuye aproximadamente de manera normal con una varianza correspondiente a la cantidad que varía el auxiliar con respecto a las variables normales (EDITAR: como lo señaló Alecos Papadopoulos, el método delta puede expresarse de manera más general, de modo que no requiere normalidad asintótica). La forma más fácil de pensar en esto es como una expansión de Taylor, donde el primer término de una función es la media, y la varianza proviene de los términos de segundo orden. En concreto, si es una función del parámetro β y b es un estimador consistente, normalmente distribuido para ese parámetro: g $g$$\beta$$b$ Dado que β es una constante y b es un estimador consistente para β , podemos decir:

$$g(b) \approx g(\beta) + \nabla g(\beta)^\prime (b - \beta)$$ $\beta$$b$$\beta$ En este caso, b es su estimación de MCO y g es el AME. Puede escribir este AME específico como: g ( b 1 , b 2 ) = b 1 + b 2  media ( x 2 ) si tomó el gradiente de esta función (recuerde, una función de loscoeficientesno de $$\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)$$ $b$$g$ $$g(b_1,b_2)=b_1+b_2 \text{ mean}(x_2)$$ ), sería: [ 1 $x_2$ y la matriz de varianza-covarianza para b podría ser: [ s 11 s 12 s 12 s 22 ] Al conectar esto a la fórmula de varianza y hacer un poco de álgebra matricial, obtendrá la misma expresión que deseaba. $$[1,\,\, \text{mean}(x_2)]^\prime$$ $b$ $$\left[ \begin{matrix} s_{11} & s_{12} \\ s_{12} & s_{22} \end{matrix}\right]$$

$g$RnumDeriv

ADENDA: En este caso específico, el Rcódigo sería:

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 

$g$