如何在R的spdep包中为空间滞后模型（lagsarlm）添加回归权重？

问题背景 #

在使用R的spdep包估计空间滞后模型（SLM）时，发现lagsarlm函数不支持weights参数，而同包的errorsarlm（空间误差模型SEM）却支持该参数。需求是基于县级别聚合数据（选举参与率+社会人口变量），按各县选民数量加权，同时考虑空间滞后效应（邻近县参与率的影响），模型形式为：
$$y=\rho Wy+X\beta+\varepsilon$$
其中$W$为空间距离权重矩阵，$\rho$为空间滞后系数，$\beta$为自变量系数，$\varepsilon$为随机误差项。

为什么lagsarlm不支持weights参数？ #

核心原因在于两种空间模型的估计逻辑差异：

• 空间误差模型（SEM）：通常通过广义最小二乘（GLS）估计，GLS天然支持加权操作——只需将权重纳入协方差矩阵的调整中，errorsarlm的实现整合了这一逻辑。
• 空间滞后模型（SLM）：属于内生模型（空间滞后项$Wy$与误差项$\varepsilon$相关），spdep中lagsarlm采用极大似然（ML）或工具变量（IV）估计，但当前版本未整合加权逻辑。加权SLM的ML估计需要修改对数似然函数，纳入权重对残差平方和的加权项，这一开发在现有lagsarlm中未完成。

实现加权空间滞后模型的替代方案 #

方案1：变量加权变换后使用lagsarlm #

通过变量变换将加权问题转化为标准SLM的估计，思路类似加权最小二乘：将因变量和自变量分别乘以权重的平方根，再用lagsarlm估计变换后的模型。

library(spdep)

# 加载示例数据
example(columbus)
listw <- nb2listw(col.gal.nb)

# 生成模拟权重（如选民数量）
myweights <- rnorm(nrow(columbus), mean = 1000, sd = 100)
# 确保权重为正（避免负数影响）
myweights <- pmax(myweights, 1)
sqrt_w <- sqrt(myweights)

# 构造加权后的数据集
columbus_w <- columbus
columbus_w$CRIME <- columbus$CRIME * sqrt_w
columbus_w$INC <- columbus$INC * sqrt_w
columbus_w$HOVAL <- columbus$HOVAL * sqrt_w

# 估计加权空间滞后模型
w_slag <- lagsarlm(CRIME ~ INC + HOVAL, data = columbus_w, listw = listw)
summary(w_slag)

方案2：自定义极大似然估计函数 #

如果需要更严格的加权ML估计，可以手动构造对数似然函数，通过optim工具完成参数优化：

# 定义加权SLM的对数似然函数（取负用于optim最小化）
loglik_weighted_slm <- function(par, y, X, W, weights) {
  rho <- par[1]
  beta <- par[-1]
  n <- length(y)
  
  # 计算残差
  e <- y - rho * as.vector(W %*% y) - X %*% beta
  # 加权残差平方和
  weighted_ss <- sum(weights * e^2)
  # 方差的MLE
  sigma2 <- weighted_ss / n
  
  # 对数似然值（取负）
  ll <- -n/2*log(2*pi) - n/2*log(sigma2) - weighted_ss/(2*sigma2)
  return(-ll)
}

# 准备数据
y <- columbus$CRIME
# 提取自变量矩阵（去除截距项）
X <- model.matrix(CRIME ~ INC + HOVAL, data = columbus)[, -1]
# 转换空间权重为矩阵形式
W <- as.matrix(listw)

# 初始参数值（rho初始为0，beta初始为普通OLS结果）
init_par <- c(0, coef(lm(CRIME ~ INC + HOVAL, data = columbus)))

# 优化求解
opt_result <- optim(init_par, loglik_weighted_slm, 
                    y = y, X = X, W = W, weights = myweights,
                    method = "BFGS")

# 提取并输出估计结果
rho_est <- opt_result$par[1]
beta_inc <- opt_result$par[2]
beta_hoval <- opt_result$par[3]
sigma2_est <- sum(myweights * (y - rho_est * as.vector(W %*% y) - X %*% c(beta_inc, beta_hoval))^2)/nrow(columbus)

cat("加权空间滞后模型手动ML估计结果：\n")
cat("空间滞后系数 rho:", round(rho_est, 4), "\n")
cat("INC系数 beta:", round(beta_inc, 4), "\n")
cat("HOVAL系数 beta:", round(beta_hoval, 4), "\n")
cat("误差项方差 sigma²:", round(sigma2_est, 4), "\n")

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