Rの関数:lagsarlm {spatialreg }

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

# 乱数シードの固定（再現性確保のため）
seed <- 20251127
set.seed(seed)

# -----------------------------------------------------------------------------
# 1. 空間構造とサンプルデータの作成
# -----------------------------------------------------------------------------

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

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

# 近傍リスト（neighbours list）の作成：クイーン基準
# 縦・横・斜めの接するセルを隣接とみなします
nb <- cell2nb(side_length, side_length, type = "queen")

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

# データ生成計算用に、空間重み行列を疎行列形式（Sparse Matrix）に変換
W_mat <- as(listw_obj, "CsparseMatrix")

# -----------------------------------------------------------------------------
# 2. 真のモデルパラメータの設定
# -----------------------------------------------------------------------------

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

# 真の回帰係数 (beta)
beta_0 <- 2.0 # 切片
beta_1 <- 1.5 # X1の係数
beta_2 <- -0.8 # X2の係数

# 真の空間ラグパラメータ (rho)
rho_true <- 0.6

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

# -----------------------------------------------------------------------------
# 3. 従属変数 y の生成 (空間ラグモデルのデータ生成プロセス)
# -----------------------------------------------------------------------------

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

I_mat <- Diagonal(n)

# (I - rho * W) の逆行列を計算
# rhoが存在することで、yの値は隣接地域のyの値に依存することになります
inv_mat <- solve(I_mat - rho_true * W_mat)

# 線形予測子 (X * beta)
X_beta <- beta_0 + beta_1 * X1 + beta_2 * X2

# 従属変数 y を生成
y_vec <- as.vector(inv_mat %*% (X_beta + epsilon))

# 解析用データフレームの作成
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))

# -----------------------------------------------------------------------------
# 6. 結果の表示と考察
# -----------------------------------------------------------------------------

cat("----------------------------------------------------------------\n")
cat("【推定結果の比較】\n")
cat("真のパラメータ:\n")
cat(sprintf(
  "  (Intercept): %.2f\n  x1:          %.2f\n  x2:          %.2f\n  rho:         %.2f\n",
  beta_0, beta_1, beta_2, rho_true
))
cat("----------------------------------------------------------------\n")

cat("\n(A) 通常の最小二乗法 (OLS) の結果:\n")
# OLSの係数表を表示（空間相関のバイアスが含まれる可能性があります）
print(summary(model_ols)$coefficients[, 1:2])

cat("\n(B) 空間ラグモデル (lagsarlm) の結果:\n")
# summaryオブジェクトの取得
sum_sar <- summary(model_sar)

# 係数表の表示
# sum_sar$Coef を使用して、推定値と標準誤差を抽出します
print(sum_sar$Coef[, 1:2])

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