在R中为SAR模型绘制带预测区间的部分依赖图

核心思路 #

由于predict.sarlm不支持直接输出置信区间，我们通过参数Bootstrap模拟来捕获参数不确定性，进而生成预测值的分布并计算95%置信区间，最终绘制带区间的部分依赖图。

步骤1：准备基础数据与模型 #

先拟合SAR模型，同时构建固定其他协变量为均值、仅目标变量取序列值的预测数据集：

library(spdep)
library(tidyverse)
library(MASS)

# 加载示例数据并拟合SAR模型
data(oldcol)
lw <- nb2listw(COL.nb, style="W")
mod <- errorsarlm(CRIME ~ INC + HOVAL, data=COL.OLD, lw)

# 提取协变量均值，用于固定其他变量
cov_means <- colMeans(COL.OLD[, c("INC", "HOVAL")])

# 生成目标变量（以INC为例）的取值序列，覆盖原数据范围
inc_seq <- seq(min(COL.OLD$INC), max(COL.OLD$INC), length.out = 50)

# 构建预测数据集：固定HOVAL为均值，INC取序列值
pred_data_inc <- tibble(
  INC = inc_seq,
  HOVAL = rep(cov_means["HOVAL"], length(inc_seq))
)

步骤2：参数Bootstrap模拟预测值 #

从模型参数的渐近正态分布中抽样，生成多组参数并计算对应预测值：

# 获取模型参数的均值和方差-协方差矩阵
param_mean <- coef(mod)
param_vcov <- vcov(mod)

# 设置Bootstrap次数（建议至少1000次，示例用200次简化计算）
n_boot <- 200

# 初始化存储预测值的矩阵
boot_preds_inc <- matrix(NA, nrow = length(inc_seq), ncol = n_boot)

# 固定随机种子保证结果可重复
set.seed(123)
for (i in 1:n_boot) {
  # 从参数的正态分布中抽样
  boot_params <- mvrnorm(1, param_mean, param_vcov)
  
  # 替换模型参数并计算预测值
  boot_mod <- mod
  boot_mod$coefficients <- boot_params
  # 若需同时考虑lambda的不确定性，可从其渐近分布抽样替换此处固定值
  boot_mod$lambda <- mod$lambda
  
  boot_preds_inc[, i] <- predict(boot_mod, newdata = pred_data_inc, listw = lw)
}

步骤3：计算置信区间并绘图 #

基于Bootstrap预测值分布，计算95%置信区间并绘制部分依赖图：

# 计算每个INC取值对应的预测均值、2.5%和97.5%分位数（置信区间）
pred_stats_inc <- pred_data_inc %>%
  mutate(
    pred_mean = rowMeans(boot_preds_inc),
    pred_lower = apply(boot_preds_inc, 1, quantile, 0.025),
    pred_upper = apply(boot_preds_inc, 1, quantile, 0.975)
  )

# 绘制带95%置信区间的部分依赖图
ggplot(pred_stats_inc, aes(x = INC, y = pred_mean)) +
  geom_line(color = "darkblue", linewidth = 1) +
  geom_ribbon(aes(ymin = pred_lower, ymax = pred_upper), fill = "lightblue", alpha = 0.3) +
  labs(x = "家庭收入（INC）", y = "犯罪率预测值", title = "INC对CRIME的部分依赖图（带95%置信区间）") +
  theme_minimal()

针对HOVAL的重复操作 #

只需将目标变量替换为HOVAL，重复上述步骤即可：

# 生成HOVAL的取值序列
hoval_seq <- seq(min(COL.OLD$HOVAL), max(COL.OLD$HOVAL), length.out = 50)
pred_data_hoval <- tibble(
  INC = rep(cov_means["INC"], length(hoval_seq)),
  HOVAL = hoval_seq
)

# 后续Bootstrap模拟、置信区间计算、绘图步骤与INC完全一致，替换pred_data即可

注意事项 #

• 若追求更高精度，建议将n_boot调整至1000次以上，计算时间会相应增加。
• 若需考虑空间自回归系数lambda的不确定性，可通过mod$se$lambda获取其标准误，从正态分布中抽样替换步骤2中固定的lambda值。
• 若需控制其他变量为中位数而非均值，修改cov_means的计算方式为apply(COL.OLD[, c("INC", "HOVAL")], 2, median)即可。

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