Rの関数:errorsarlm {spatialreg}

Rの関数から errorsarlm {spatialreg} を確認します。

本ポストはこちらの続きです。

const typesetMath = (el) => { if (window.MathJax) { // MathJax Typeset window.MathJax.typeset(); } else if (window.katex...

関数 `errorsarlm` とは

errorsarlm（Spatial Error Model: SEM）は、空間誤差モデルを最尤法を用いて推定するための関数です。

空間誤差モデルは、従属変数 \(y\) ではなく、誤差項 \(u\)に空間的な自己相関が存在すると仮定するモデルです。

観測できない要因が空間的に相関している場合に利用されます。

数式では以下のように表現されます。

\[\begin{eqnarray}
y &=& X \beta + u \\
u &=& \lambda W u + \epsilon
\end{eqnarray}\]

ここで、\(u\)は空間的自己相関を持つ誤差項、\(\lambda\)は空間誤差パラメータ、\(\epsilon\)は独立した誤差項（ホワイトノイズ）です。

関数 `errorsarlm` の引数

library(spatialreg)
args(errorsarlm)

function (formula, data = list(), listw, na.action, weights = NULL, 
    Durbin, etype, method = "eigen", quiet = NULL, zero.policy = NULL, 
    interval = NULL, tol.solve = .Machine$double.eps, trs = NULL, 
    control = list()) 
NULL

formula
- モデルの構造を指定する式です（例: y ~ x1 + x2）。通常の回帰モデルと同様に記述します。
data
- formulaに含まれる変数が格納されているデータフレームです。デフォルトは空のリストですが、分析対象のデータを指定します。
listw
- 空間重み行列を格納したlistwクラスのオブジェクトです。誤差項の空間的波及構造 \(W\) を定義します。
na.action
- 欠損値（NA）の取り扱いを指定する関数です。欠損があるデータの行を除外するかどうかを決定します。
weights
- 重み付き最小二乗法（WLS）を行う際の観測値ごとの重みを指定する数値ベクトルです。分散が均一でない場合などに使用されます。デフォルトはNULL（すべての重みが1）です。
Durbin
- 空間ダービン誤差モデル（Spatial Durbin Error Model: SDEM）を推定するかを指定します。
- TRUEの場合、説明変数の空間ラグ（\(WX\)）がモデルに追加されます（\(y = X\beta + WX\gamma + u\)）。
- デフォルトはetypeが"error"の場合はFALSEとなります。
etype
- 誤差モデルのタイプを指定します。
  - "error": 通常の空間誤差モデル（SEM）。
  - "emixed": 空間ダービン誤差モデル（SDEM）。Durbin引数が指定された場合、自動的にこちらが選択されます。
method
- ヤコビアンの計算方法を指定します。デフォルトは"eigen"（固有値を用いる方法）です。"Chebyshev"や"LU"などの近似・高速化手法が選択可能です。
quiet
- 実行時のメッセージ出力を抑制するかどうかを指定する論理値です。
zero.policy
- 隣接地点を持たない（孤立した）観測地点を許容するかを指定します。TRUEの場合、計算を続行します。
interval
- 空間誤差パラメータ \(\lambda\) の探索範囲を指定する数値ベクトルです。
tol.solve
- 行列演算（逆行列計算等）における許容誤差（tolerance）です。数値的な安定性を確保するために使用されます。
trs
- 重み行列 \(W\) のトレース（対角成分の和）の近似値を指定します。近似手法を用いる際の計算コスト削減のために利用されます。
control
- 最適化アルゴリズムや数値計算の詳細設定をリスト形式で渡します。ヘッセ行列の計算有無や、収束判定基準などを微調整します。

シミュレーション用サンプルデータとコード

errorsarlmの機能を確認するため、以下の手順でサンプルデータを作成し、シミュレーションを行います。

空間構造の生成
- \(20 \times 20\) のグリッド空間（計400地点）を作成します。
真のパラメータ設定
- 回帰係数 \(\beta\) と、誤差項の空間依存パラメータ \(\lambda\) を設定します。
データ生成:
- まず、独立した誤差項 \(\epsilon\) を生成します。
- 次に、空間的自己相関を持つ誤差項 \(u = (I - \lambda W)^{-1} \epsilon\) を生成します。
- 最後に、従属変数 \(y = X \beta + u\) を計算します。
比較検証
- 通常の線形回帰（OLS）と errorsarlm の結果を比較します。

以下はRコードです。

空間構造とサンプルデータの作成

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

# 乱数シードの固定
seed <- 20251129
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 = 3)
X2 <- runif(n, min = -5, max = 5)

# 真の回帰係数 (beta)
beta_0 <- 5.0 # 切片
beta_1 <- 2.5 # X1の係数
beta_2 <- -1.2 # X2の係数

# 真の空間誤差パラメータ (lambda)
# 誤差項のみに空間的な相関を持たせます
lambda_true <- 0.7

# 独立した誤差項 (epsilon): ホワイトノイズ
epsilon <- rnorm(n, mean = 0, sd = 2)

従属変数 y の生成 (空間誤差モデルのデータ生成プロセス)

空間誤差モデル（SEM）では誤差項 \(u = (I - \lambda W)^{-1}\epsilon\) を逆行列を用いて生成します。

\(y\) は \(X\beta\) （構造部分）にこの空間相関のある誤差 \(u\) を単純に加算して作成されます。

このデータ生成構造が errorsarlm の前提と一致します。

# モデル式: y = X * beta + u
# 誤差構造: u = lambda * W * u + epsilon
# 変形: (I - lambda * W) * u = epsilon
# 解: u = (I - lambda * W)^(-1) * epsilon

I_mat <- Diagonal(n)

# 空間的自己相関を持つ誤差項 u を生成
inv_mat <- solve(I_mat - lambda_true * W_mat)
u_vec <- as.vector(inv_mat %*% epsilon)

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

# 従属変数 y を生成 (構造部分 + 空間相関のある誤差)
y_vec <- X_beta + u_vec

# 解析用データフレームの作成
sim_data <- data.frame(
  id = 1:n,
  x_coord = coords$x,
  y_coord = coords$y,
  y = y_vec,
  x1 = X1,
  x2 = X2,
  u_true = u_vec # 可視化確認用の真の誤差項
)

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

--- 生成したデータの一部を確認 ---
'data.frame':   400 obs. of  7 variables:
 $ id     : int  1 2 3 4 5 6 7 8 9 10 ...
 $ x_coord: int  1 2 3 4 5 6 7 8 9 10 ...
 $ y_coord: int  1 1 1 1 1 1 1 1 1 1 ...
 $ y      : num  29.9 38.6 29 25.9 33.1 ...
 $ x1     : num  10.13 11.41 10.39 5.75 8.93 ...
 $ x2     : num  1.87 -4.81 1.3 -4.81 -4.14 ...
 $ u_true : num  1.783 -0.72 -0.41 0.781 0.784 ...
NULL

データの可視化 (誤差項の空間分布)

# 真の誤差項 u の空間分布プロット
# 説明変数で説明できない部分に、クラスター（空間相関）があることを確認します
p1 <- ggplot(sim_data, aes(x = x_coord, y = y_coord, fill = u_true)) +
  geom_tile() +
  scale_fill_viridis_c(option = "mako") +
  labs(
    title = "シミュレーションデータ: 真の誤差項 u の空間分布",
    subtitle = paste("真の空間誤差パラメータ lambda =", lambda_true),
    x = "X 座標", y = "Y 座標", fill = "誤差値"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(face = "bold"))

print(p1)

視覚的特徴と色の意味

横軸（X座標）と縦軸（Y座標）は \(20 \times 20\) のグリッド上の位置を表します。
右側の凡例が示す通り、値が低い（マイナス）領域は濃い紫色・黒色で、値が高い（プラス）領域は明るい緑色・白色で表現されています。

空間的自己相関の現れ（クラスタリング）

本シミュレーションでは、空間誤差パラメータ \(\lambda\)（lambda）を 0.7 という比較的高い値に設定しています。Figure 1 を確認しますと、色がランダムな砂嵐（ホワイトノイズ）のようになっておらず、似た色が塊（クラスター）を作っています。
ある地点の誤差項が高い値であれば、その隣接地点も高い値を取りやすいという「正の空間的自己相関」が確認できます。こちらは、数式 \(u = (I - \lambda W)^{-1}\epsilon\) によって、隣接地域間でショックが波及・共有された結果です。

モデルの推定

引数 listw に空間重み行列を指定することで、関数内部で誤差項の共分散構造 \(\Omega = \sigma^2 (I - \lambda W)^{-1}(I - \lambda W')^{-1}\) を考慮した最尤推定が行われます。

# (A) 通常の線形回帰モデル (OLS)
# 誤差項の空間相関を無視します
model_ols <- lm(y ~ x1 + x2, data = sim_data)

# (B) 空間誤差モデル (errorsarlm)
# 誤差項の空間構造をモデル化します
model_sem <- errorsarlm(y ~ x1 + x2, data = sim_data, listw = listw_obj)

cat("--- 推定された線形回帰モデル ---\n")
summary(model_ols)

cat("--- 推定された空間誤差モデル ---\n")
summary(model_sem)

--- 推定された線形回帰モデル ---

Call:
lm(formula = y ~ x1 + x2, data = sim_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.9809 -1.5897 -0.0345  1.6589  6.1917 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.80697    0.41081   11.70   <2e-16 ***
x1           2.53370    0.03934   64.40   <2e-16 ***
x2          -1.19124    0.04047  -29.43   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.314 on 397 degrees of freedom
Multiple R-squared:  0.927, Adjusted R-squared:  0.9266 
F-statistic:  2521 on 2 and 397 DF,  p-value: < 2.2e-16

--- 推定された空間誤差モデル ---

Call:errorsarlm(formula = y ~ x1 + x2, data = sim_data, listw = listw_obj)

Residuals:
      Min        1Q    Median        3Q       Max 
-5.780925 -1.355678 -0.021658  1.446636  5.862258 

Type: error 
Coefficients: (asymptotic standard errors) 
             Estimate Std. Error z value  Pr(>|z|)
(Intercept)  5.154126   0.428349  12.033 < 2.2e-16
x1           2.495894   0.032533  76.719 < 2.2e-16
x2          -1.206412   0.033388 -36.133 < 2.2e-16

Lambda: 0.64816, LR test value: 97.152, p-value: < 2.22e-16
Asymptotic standard error: 0.054409
    z-value: 11.913, p-value: < 2.22e-16
Wald statistic: 141.91, p-value: < 2.22e-16

Log likelihood: -853.1296 for error model
ML residual variance (sigma squared): 3.8627, (sigma: 1.9654)
Number of observations: 400 
Number of parameters estimated: 5 
AIC: 1716.3, (AIC for lm: 1811.4)

モデルの比較

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

cat("\n(A) 通常の最小二乗法 (OLS) の結果:\n")
# 係数と標準誤差を表示
print(summary(model_ols)$coefficients[, 1:2])

cat("\n(B) 空間誤差モデル (errorsarlm) の結果:\n")
# summaryオブジェクトの取得
sum_sem <- summary(model_sem)

# 推定された係数表 (Coefを使用)
print(sum_sem$Coef[, 1:2])

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

【推定結果の比較】
真のパラメータ:
  (Intercept): 5.00
  x1:          2.50
  x2:          -1.20
  lambda:      0.70
----------------------------------------------------------------

(A) 通常の最小二乗法 (OLS) の結果:
             Estimate Std. Error
(Intercept)  4.806974 0.41081451
x1           2.533699 0.03934274
x2          -1.191240 0.04047418

(B) 空間誤差モデル (errorsarlm) の結果:
             Estimate Std. Error
(Intercept)  5.154126 0.42834858
x1           2.495894 0.03253308
x2          -1.206412 0.03338813

推定された空間誤差パラメータ lambda: 0.6482 (標準誤差: 0.0544)

回帰係数（Estimate）の比較

まず、OLSとSEMの係数推定値を比較します。

(A) OLS: x1 = 2.53, x2 = -1.19
(B) SEM: x1 = 2.50, x2 = -1.21

空間誤差モデル（SEM）の設定下では、OLS推定量も統計的に「不偏かつ一致推定量」となります。

すなわち、誤差項に空間相関があっても、OLSによる係数の推定値自体に大きなバイアス（歪み）は生じません。

結果として、OLSもSEMも真の値（2.50, -1.20）に近い値を推定しています。

標準誤差（Std. Error）と効率性

次に、推定の精度を示す「標準誤差」を確認します。

x1の標準誤差: OLS (0.039) > SEM (0.033)
x2の標準誤差: OLS (0.040) > SEM (0.033)

OLSは「誤差項が互いに独立である」と仮定して計算を行いますが、実際には空間相関があるため、情報を効率的に使えません。

対して、errorsarlm は空間重み行列 \(W\) を用いて誤差の相関構造をモデルに組み込み、一般化最小二乗法（GLS）に近い処理を行います。

その結果、SEMのほうが標準誤差が小さくなっています。

標準誤差が小さいということは、「推定値のばらつきが小さく、より信頼性が高い（効率的である）」ことを意味します。

空間誤差パラメータ (lambda) の特定

推定値: 0.6482
真の値: 0.70

推定された \(\lambda\) は 0.6482 であり、標準誤差（0.0544）を考慮した95%信頼区間の中に真の値 0.70 が含まれています。

下限: \(0.6482 - (0.0544 \times 2) = 0.6482 - 0.1088 = \mathbf{0.5394}\)
上限: \(0.6482 + (0.0544 \times 2) = 0.6482 + 0.1088 = \mathbf{0.7570}\)

OLSではこの「観測されない空間的なパターン（誤差の相関）」を数値化することはできませんが、errorsarlm を用いることで、背景にある空間的な構造を定量的に特定することが可能になります。

以上です。

関数 errorsarlm とは

関数 errorsarlm の引数