Rの関数:lagsarlm(type = “mixed”) {spatialreg }

# 必要なパッケージの読み込み
library(spdep)
library(spatialreg)
library(Matrix)
library(ggplot2)

# 乱数シードの固定
seed <- 20251201
set.seed(seed)

# グリッドサイズの設定（20x20 = 400地点）
side_length <- 20
n <- side_length * side_length

# 座標データの作成
coords <- expand.grid(x = 1:side_length, y = 1:side_length)

# 近傍リストの作成（クイーン基準）
nb <- cell2nb(side_length, side_length, type = "queen")

# 空間重み行列（listwオブジェクト）の作成：行規格化
listw_obj <- nb2listw(nb, style = "W")
W_mat <- as(listw_obj, "CsparseMatrix")

# 説明変数 X1, X2 の生成
X1 <- rnorm(n, mean = 10, sd = 2)
X2 <- runif(n, min = 0, max = 5)

# 説明変数の空間ラグ (WX) を計算
# 隣接地域の説明変数の平均値が、自地域に影響を与えると仮定します
WX1 <- as.vector(W_mat %*% X1)
WX2 <- as.vector(W_mat %*% X2)

# 回帰係数 (beta): 自地域の説明変数が及ぼす直接的な影響
beta_0 <- 3.0 # 切片
beta_1 <- 1.5 # X1の係数
beta_2 <- -2.0 # X2の係数

# 空間ラグ係数 (gamma): 隣接地域の説明変数が及ぼす影響
gamma_1 <- 0.8 # WX1の係数 (正のスピルオーバー)
gamma_2 <- -1.5 # WX2の係数 (負のスピルオーバー)

# 空間自己回帰パラメータ (rho): 隣接地域の従属変数が及ぼす影響
rho_true <- 0.5

# 誤差項 (epsilon)
epsilon <- rnorm(n, mean = 0, sd = 1)

# モデル式: y = rho * Wy + X * beta + WX * gamma + epsilon
# 変形: (I - rho * W) * y = X * beta + WX * gamma + epsilon
# 解: y = (I - rho * W)^(-1) * (X * beta + WX * gamma + epsilon)

# 右辺の構造部分 (Trend + Lagged X + Error)
# ここで WX1, WX2 の項を含めて SDM のデータ生成プロセスとします
X_terms <- beta_0 + beta_1 * X1 + beta_2 * X2
WX_terms <- gamma_1 * WX1 + gamma_2 * WX2
total_disturbance <- X_terms + WX_terms + epsilon

# 逆行列による y の算出
I_mat <- Diagonal(n)
inv_mat <- solve(I_mat - rho_true * W_mat)
y_vec <- as.vector(inv_mat %*% total_disturbance)

# 解析用データフレームの作成
sim_data <- data.frame(
  id = 1:n,
  x_coord = coords$x,
  y_coord = coords$y,
  y = y_vec,
  x1 = X1,
  x2 = X2
)

# 生成したサンプルデータの一部を確認
cat("--- 生成したサンプルデータの一部を確認 ---\n")
print(str(sim_data))

# lagsarlm 関数で type = "mixed" を指定
# こちらも Durbin = TRUE と指定することと同義です
model_sdm <- lagsarlm(y ~ x1 + x2, data = sim_data, listw = listw_obj, type = "mixed")

summary(model_sdm)

cat("----------------------------------------------------------------\n")
cat("【推定結果の比較: 空間ダービンモデル (SDM)】\n")
cat("----------------------------------------------------------------\n")
cat("真のパラメータ設定:\n")
cat(sprintf("  (Intercept): %.2f\n", beta_0))
cat(sprintf("  x1 (beta):   %.2f\n", beta_1))
cat(sprintf("  x2 (beta):   %.2f\n", beta_2))
cat(sprintf("  lag.x1 (gamma): %.2f\n", gamma_1))
cat(sprintf("  lag.x2 (gamma): %.2f\n", gamma_2))
cat(sprintf("  rho:         %.2f\n", rho_true))
cat("----------------------------------------------------------------\n")

cat("\n(B) lagsarlm(type = 'mixed') の推定結果:\n")
# summaryオブジェクトの取得
sum_sdm <- summary(model_sdm)

# 推定された係数表 (Coef)
# ここには X の係数だけでなく、lag.X (WX) の係数も表示されます
print(sum_sdm$Coef[, 1:2])

cat(sprintf(
  "\n推定された空間自己回帰パラメータ rho: %.4f (標準誤差: %.4f)\n",
  model_sdm$rho, model_sdm$rho.se
))