Bagaimana cara menggunakan metode delta untuk kesalahan standar efek marginal?

Saya tertarik untuk lebih memahami metode delta untuk mendekati kesalahan standar dari efek marginal rata-rata dari model regresi yang mencakup istilah interaksi. Saya telah melihat pertanyaan terkait di bawah metode delta tetapi tidak ada yang memberikan apa yang saya cari.

Pertimbangkan contoh data berikut ini sebagai contoh yang memotivasi:

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

Saya tertarik dengan efek marginal rata-rata (AME) dari x1dan x2. Untuk menghitung ini, saya cukup melakukan hal berikut:

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

Tetapi bagaimana saya menggunakan metode delta untuk menghitung kesalahan standar AME ini?

Saya dapat menghitung SE untuk interaksi khusus ini dengan tangan:

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

Tapi saya tidak mengerti cara menggunakan metode delta.

Idealnya, saya mencari beberapa panduan tentang cara berpikir tentang (dan kode) metode delta untuk AME dari setiap model regresi sewenang-wenang. Misalnya, pertanyaan ini memberikan rumus untuk SE untuk efek interaksi tertentu dan dokumen ini dari Matt Golder menyediakan rumus untuk berbagai model interaktif, tetapi saya ingin lebih memahami prosedur umum untuk menghitung UK dari AMEs daripada rumus untuk SE dari AME tertentu.

Metode delta hanya mengatakan bahwa jika Anda dapat mewakili variabel bantu Anda dapat mewakili sebagai fungsi dari variabel acak yang terdistribusi normal, bahwa variabel tambahan terdistribusi secara normal dengan varians yang sesuai dengan seberapa banyak variasi pembantu terkait dengan variabel normal (EDIT: seperti yang ditunjukkan oleh Alecos Papadopoulos, metode delta dapat dinyatakan lebih umum sehingga tidak memerlukan normalitas asimptotik). Cara termudah untuk memikirkan ini adalah sebagai ekspansi Taylor, di mana istilah pertama suatu fungsi adalah rata-rata, dan varians berasal dari istilah urutan kedua. Secara khusus, jika adalah fungsi dari parameter dan adalah penduga yang konsisten dan terdistribusi normal untuk parameter itu: $g$$\beta$$b$

$$g(b) \approx g(\beta) + \nabla g(\beta)^\prime (b - \beta)$$ Karena adalah konstanta, dan adalah estimator yang konsisten untuk , kita dapat mengatakan: Dalam hal ini, adalah taksiran OLS Anda, dan adalah AME. Anda dapat menulis AME spesifik ini sebagai: jika Anda mengambil gradien dari fungsi ini (ingat, fungsi dari koefisien bukan dari ), itu akan menjadi : $\beta$$b$$\beta$ b g g ( b 1 , b 2 ) = b 1 + b 2  rata-rata ( x 2 ) x 2 [ 1 $$\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$ b [ s 11 s 12 s 12 s 22 $$[1,\,\, \text{mean}(x_2)]^\prime$$ dan matriks varians-kovariansi untuk mungkin: Memasukkan ini ke rumus varians dan melakukan beberapa aljabar matriks memberi Anda ekspresi yang sama yang Anda inginkan.$b$ $$\left[ \begin{matrix} s_{11} & s_{12} \\ s_{12} & s_{22} \end{matrix}\right]$$

Secara umum jika Anda ingin melakukan ini, Anda dapat secara eksplisit kode apapun Anda ingin menjadi sebagai fungsi dari semua koefisien Anda dan kemudian menggunakan untuk mengambil gradien numerik (jika tidak, anda harus menggunakan komputer aljabar) fungsi sehubungan dengan parameter Anda, pada parameter yang Anda perkirakan. Maka Anda cukup mengambil varians-kovarians matriks dan gradien numerik ini dan hubungkan ke formula dan voila! Metode delta.$g$RnumDeriv

ADDENDUM: Dalam kasus khusus ini Rkodenya adalah:

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 

Perhatikan bahwa akan selalu lebih disukai untuk mendapatkan gradien yang tepat daripada gradien numerik untuk masalah ini, karena gradien yang tepat akan memiliki lebih sedikit kesalahan komputasi. Fakta bahwa linier menghilangkan masalah ini, dan untuk fungsi yang lebih rumit, gradien yang tepat mungkin tidak selalu tersedia.$g$