Jak stosować metodę delta do standardowych błędów efektów marginalnych?

Interesuje mnie lepsze zrozumienie metody delta do aproksymacji standardowych błędów średnich efektów krańcowych modelu regresji, który obejmuje termin interakcji. Patrzyłem na powiązane pytania metodą delta, ale żadne nie dostarczyło tego, czego szukam.

Rozważ następujące przykładowe dane jako motywujący przykład:

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

Interesują mnie średnie efekty krańcowe (AME) x1i x2. Aby je obliczyć, po prostu wykonuję następujące czynności:

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

Ale jak użyć metody delta do obliczenia standardowych błędów tych AME?

Mogę ręcznie obliczyć SE dla tej konkretnej interakcji:

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

Ale nie rozumiem, jak używać metody delta.

Idealnie szukam wskazówek, jak myśleć (i kodować) metodę delta dla AME dowolnego modelu regresji arbitralnej. Na przykład to pytanie zawiera formułę dla SE dla określonego efektu interakcji, a ten dokument Matta Goldera zawiera formuły dla różnych modeli interaktywnych, ale chcę lepiej zrozumieć ogólną procedurę obliczania SE dla AME, a nie formułę dla SE dowolnego konkretnego lekarza orzecznika medycyny lotniczej.

Metoda delta mówi po prostu, że jeśli możesz reprezentować zmienną pomocniczą, możesz reprezentować ją jako funkcję normalnie rozmieszczonych zmiennych losowych, ta zmienna pomocnicza jest w przybliżeniu normalnie rozmieszczona z wariancją odpowiadającą, jak bardzo zmienna pomocnicza zmienia się względem normalnych zmiennych (EDYCJA: jak wskazał Alecos Papadopoulos, metodę delta można określić bardziej ogólnie, że nie wymaga ona asymptotycznej normalności). Najłatwiejszym sposobem na myślenie o tym jest rozwinięcie Taylora, w którym pierwszy składnik funkcji jest średnią, a wariancja pochodzi od warunków drugiego rzędu. W szczególności, jeśli jest funkcją parametru a jest spójnym, normalnie rozłożonym estymatorem dla tego parametru: β b g ( b ) ≈ g ( β ) + ∇ g ( β ) ′ ( b - β ) β b β $g$$\beta$$b$

$$g(b) \approx g(\beta) + \nabla g(\beta)^\prime (b - \beta)$$ Ponieważ jest stałą, a$\beta$$b$ jest spójnym estymatorem dla , możemy zatem powiedzieć: W tym przypadku jest oszacowaniem OLS, a AME. Możesz zapisać ten konkretny AME jako: jeśli wziąłeś gradient tej funkcji (pamiętaj, funkcja współczynników nie ), to będzie : $\beta$ $$\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)$$ $x_2$ $$[1,\,\, \text{mean}(x_2)]^\prime$$ a macierz wariancji-kowariancji dla może wyglądać następująco: Podłączenie Formuła wariancji i wykonanie algebry macierzowej daje takie samo wyrażenie, jak chcesz.$b$ $$\left[ \begin{matrix} s_{11} & s_{12} \\ s_{12} & s_{22} \end{matrix}\right]$$

W ogóle, jeśli chcesz to zrobić, można jawnie kod cokolwiek chcesz się jako funkcja wszystkich współczynników, a następnie użyć do podjęcia gradient liczbową (w przeciwnym razie trzeba by wykorzystanie algebry komputera) funkcji w odniesieniu do twoje parametry, według parametrów, które oszacowałeś. Następnie weź po prostu macierz wariancji-kowariancji i ten numeryczny gradient i podłącz ją do wzoru i voila! Metoda delta.$g$RnumDeriv

DODATEK: W tym konkretnym przypadku Rkod będzie:

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 

Zauważ, że zawsze lepiej będzie uzyskać dokładny gradient zamiast gradientu numerycznego dla tego problemu, ponieważ dokładny gradient będzie miał mniej błędu obliczeniowego. Fakt, że jest liniowy, eliminuje ten problem, a dla bardziej skomplikowanych funkcji dokładny gradient może nie zawsze być dostępny.$g$