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

Question

Estoy interesado en comprender mejor el método delta para aproximar los errores estándar de los efectos marginales medios de un modelo de regresión que incluye un término de interacción. He mirado preguntas relacionadas en método delta pero ninguno me ha proporcionado exactamente lo que busco.

Considere los siguientes datos como 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)

Me interesan los efectos marginales medios (EMA) de x1 et x2 . Para calcularlos, 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 utilizo el método delta para calcular los errores estándar de estas EMA?

Puedo calcular a mano el SE de esta interacción concreta:

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 utilizar 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 arbitrario. Por ejemplo, esta pregunta proporciona una fórmula para el SE de un efecto de interacción concreto y este documento de Matt Golder proporcionan fórmulas para una variedad de modelos interactivos, pero quiero entender mejor el procedimiento general para calcular los SE de los AME en lugar de la fórmula para el SE de cualquier AME en particular.

Answer 1

El método delta simplemente dice que si puedes representar una variable auxiliar que puedes representar como una función de variables aleatorias distribuidas normalmente, esa variable auxiliar está aproximadamente distribuida normalmente con la varianza correspondiente a cuánto varía la auxiliar con respecto a las variables normales. (EDIT: como ha señalado Alecos Papadopoulos, el método delta puede enunciarse de forma más general, de modo que no requiere normalidad asintótica).

La forma más sencilla de verlo es como una expansión de Taylor, en la que el primer término de una función es la media y la varianza procede de los términos de segundo orden. En concreto, si $g$ es función del parámetro $\beta$ et $b$ es un estimador consistente y normalmente distribuido para ese parámetro: $$ g(b) \approx g(\beta) + \nabla g(\beta)^\prime (b - \beta) $$ Desde $\beta$ es una constante, y $b$ es un estimador coherente de $\beta$ entonces podemos decir: $$ \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) $$ En este caso, $b$ es su estimación OLS, y $g$ es la AME. Puede escribir este AME específico como: $$ g(b_1,b_2)=b_1+b_2 \text{ mean}(x_2) $$ si tomara el gradiente de esta función (recuerde, una función de la coeficientes no de $x_2$ ), sería: $$ [1,\,\, \text{mean}(x_2)]^\prime $$ y la matriz de varianza-covarianza para $b$ podría ser: $$ \left[ \begin{matrix} s_{11} & s_{12} \\ s_{12} & s_{22} \end{matrix}\right] $$ Si lo introducimos en la fórmula de la varianza y hacemos un poco de álgebra matricial, obtendremos la misma expresión que queríamos.

En general, si quieres hacer esto, puedes codificar explícitamente lo que sea $g$ que quieras en R como una función de todos sus coeficientes y luego utilizar numDeriv para tomar el gradiente numérico (de lo contrario tendrías que usar álgebra computacional) de la función con respecto a tus parámetros, en los parámetros que estimaste. A continuación, basta con tomar la matriz de varianza-covarianza y este gradiente numérico, introducirlo en la fórmula y ¡voilá! Método Delta.

ADENDA: En este caso concreto el R có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 

Tenga en cuenta que siempre será preferible obtener el gradiente exacto en lugar del gradiente numérico para este problema, ya que el gradiente exacto tendrá menos error de cálculo. El hecho de que $g$ es lineal elimina este problema, y para funciones más complicadas el gradiente exacto puede no estar siempre disponible.