空间面板Durbin模型直接、间接与总效应估算问题咨询

我需要估算空间面板Durbin模型的直接效应、间接效应与总效应，目前没有能直接完成该估算的R包，因此手动生成了自变量的空间滞后项，代码如下：

# 生成空间滞后项
dataS$JC_2dg_s <- slag(dataS$JC_2dg, listw = cont.listw)
dataS$JD_2dg_s <- slag(dataS$JD_2dg, listw = cont.listw)

# 设定模型公式并估计
eq0_1 <- TFR_2  ~  JC_2dg  + JD_2dg  + JC_2dg_s + JD_2dg_s + 
  FLFP2030  + FLFP2030_s + logGDP_r  + logGDP_r_s + expat_ratio  + expat_ratio_s + 
  hp1000  + hp_s 
model1 <- spml(eq0_1, data=dataS, listw=cont.listw, model="within", spatial.error="none",lag=T, effect="twoway", rel.tol=2e-40)
summary(model1)

# 手动尝试计算效应
W.c <- as(as_dgRMatrix_listw(cont.listw), "CsparseMatrix")
beta_jc <- model1$coefficients[2] 
theta_jc <- model1$coefficients[4] 
beta_jd <- model1$coefficients[3]  
theta_jd <- model1$coefficients[5]

# 计算rho*W
rhoW <- model1$spat.coef * W.c    
# 单位矩阵
I <- Matrix::Diagonal(nrow(W.c))
# 计算(I - rho*W)^(-1)
inverse_matrix <- solve(I - rhoW)
# 计算beta和theta的直接影响矩阵
mateff_jc <- inverse_matrix *(beta_jc) + inverse_matrix %*% (W.c *theta_jc)
mateff_jd <- inverse_matrix *(beta_jd) + inverse_matrix %*% (W.c *theta_jd)

# 计算平均直接效应及标准误（此处存疑）
ade_jc <- mean(diag(mateff_jc))
se_ade_jc <- sd(diag(mateff_jc)) 
t_directjc <- ade_jc/se_ade_jc
p_directjc <- 2 * (1 - pt(abs(t_directjc), df = nrow(W.c) - 1))
ade_jd <- mean(diag(mateff_jd))
se_ade_jd <- sd(diag(mateff_jd)) 
t_directjd <- ade_jd/se_ade_jd
p_directjd <- 2 * (1 - pt(abs(t_directjd), df = nrow(W.c) - 1))

我查阅文献后，不确定上述标准误的计算是否正确——当前用对角线元素的标准差来估计标准误的方法忽略了模型系数估计的方差-协方差结构，这会导致标准误估计偏误。

正确的效应估算与标准误计算方法 #

1. 效应的矩阵表达式
   对于SDM模型，自变量$X_k$的总效应矩阵为：
   $$
   \text{TE}_k = (I - \rho W)^{-1}(\beta_k I + \theta_k W)
   $$

• 直接效应（Direct Effect, DE）：$\text{TE}_k$对角线元素的平均值
• 间接效应（Indirect Effect, IE）：$\text{TE}_k$所有元素的平均值减去直接效应
• 总效应（Total Effect, TE）：$\text{TE}_k$所有元素的平均值

2. 标准误的计算：Delta方法
   效应是模型系数（$\beta_k, \theta_k, \rho$）的非线性函数，需用Delta方法基于模型的方差-协方差矩阵计算标准误：

• 首先提取模型的方差-协方差矩阵
• 定义效应关于系数的梯度函数，结合方差-协方差矩阵计算效应的方差，再开平方得到标准误

以下是针对JC_2dg的修正代码示例：

# 提取模型系数及方差-协方差矩阵
coefs <- coef(model1)
vcov_mat <- vcov(model1)
rho <- model1$spat.coef
n <- nrow(W.c)

# 安装并加载numDeriv包用于计算梯度
# install.packages("numDeriv")
library(numDeriv)

# 定义计算平均直接效应的函数
calc_ade <- function(beta, theta, rho) {
  inv_mat <- solve(Matrix::Diagonal(n) - rho * W.c)
  te_mat <- inv_mat * beta + inv_mat %*% (W.c * theta)
  mean(diag(te_mat))
}

# 计算直接效应的梯度
grad_ade <- grad(func = function(par) calc_ade(par[1], par[2], par[3]),
                 x = c(coefs["JC_2dg"], coefs["JC_2dg_s"], rho))

# 计算直接效应的标准误
se_ade_jc_correct <- sqrt(t(grad_ade) %*% vcov_mat[c("JC_2dg", "JC_2dg_s", "rho"), 
                                                  c("JC_2dg", "JC_2dg_s", "rho")] %*% grad_ade)

# 计算间接效应的函数与标准误
calc_ie <- function(beta, theta, rho) {
  inv_mat <- solve(Matrix::Diagonal(n) - rho * W.c)
  te_mat <- inv_mat * beta + inv_mat %*% (W.c * theta)
  mean(te_mat) - mean(diag(te_mat))
}

grad_ie <- grad(func = function(par) calc_ie(par[1], par[2], par[3]),
                x = c(coefs["JC_2dg"], coefs["JC_2dg_s"], rho))
se_ie_jc_correct <- sqrt(t(grad_ie) %*% vcov_mat[c("JC_2dg", "JC_2dg_s", "rho"), 
                                                c("JC_2dg", "JC_2dg_s", "rho")] %*% grad_ie)

# 计算总效应的函数与标准误
calc_te <- function(beta, theta, rho) {
  inv_mat <- solve(Matrix::Diagonal(n) - rho * W.c)
  te_mat <- inv_mat * beta + inv_mat %*% (W.c * theta)
  mean(te_mat)
}

grad_te <- grad(func = function(par) calc_te(par[1], par[2], par[3]),
                x = c(coefs["JC_2dg"], coefs["JC_2dg_s"], rho))
se_te_jc_correct <- sqrt(t(grad_te) %*% vcov_mat[c("JC_2dg", "JC_2dg_s", "rho"), 
                                                c("JC_2dg", "JC_2dg_s", "rho")] %*% grad_te)

关键说明 #

• 原代码中用sd(diag(mateff_jc))计算标准误是错误的，它仅反映直接效应矩阵对角线元素的离散程度，未考虑系数估计的不确定性
• Delta方法通过计算效应函数对系数的梯度，结合系数的方差-协方差矩阵，能准确估计效应的标准误
• 对于其他自变量（如JD_2dg），只需替换对应的系数名称重复上述步骤即可

内容的提问来源于stack exchange，提问作者Chen Luo