東京大学・株式会社Nospareの菅澤です．
今回は空間データに対する回帰分析として最も有名な地理的加重回帰(GWR; Geographically Weighted Regression)について紹介します．

空間データと空間異質性

緯度・経度などの位置情報が付随したデータは一般的に空間データと呼ばれます．このようなデータは不動産価格，犯罪の発生件数，地域別の選挙投票率など，様々な場面で登場します．

具体的な問題設定として，$y_i$ $(i=1,\ldots,n)$を興味のある被説明変数とし，それに対する説明変数を$x_i$とします．(複数種類の説明変数がある場合，$x_i$はベクトルになります．)
また，各データに付随する位置情報を$s_i$とします．

通常の線形回帰分析では以下のようなモデルを想定します．

y_i=x_i^t\beta + \varepsilon_i,  \ \ \ \ i=1,\ldots,n.

このモデルの特徴は**「説明変数が被説明変数に与える効果が地点に依らず均一」**という点です．

これはどういうことを表しているのか，例えば不動産価格の例で考えてみます．

説明変数として家の敷地面積(平米)が入っていた場合，回帰係数が地点に依らず一定であることは「都心か郊外に依らず敷地面積が1平米増加したときの不動産価格の上昇度合いが変わらない」ということになります．これは都心と郊外での1平米あたりの重要性の違いを考慮しておらず，多くの場合に不適当な分析になってしまう可能性があります．

このような問題点を解決するために登場する考え方が空間異質性・空間非定常性です．端的に言うと「回帰係数が空間的に異なることを許す」考え方です．このような考えで回帰モデルを推定する有名な方法が，今回紹介する地理的加重回帰です(空間的変化係数モデルとも呼ばれています)．

地理的加重回帰(GWR)

前述のように，回帰係数$\beta$が地点ごとに異なるような以下のモデルを考えます．

y_i=x_i^t\beta_i + \varepsilon_i, \ \ \ \ \ \ i=1,\ldots,n.

$\beta_i$は$\beta(s_i)$のように位置情報の関数として表現する場合もありますが，今回は$\beta_i$という簡単な表記を用います．

各$\beta_i$を推定するには，以下のような重み付き最小二乗法を用います．

\widehat{\beta}_i={\rm argmin}_{\beta}\sum_{j=1}^n w_{ij}(y_j-x_j^t\beta)^2.  \tag{1}

ここで$w_{ij}$は$i$番目の地点を基準に$j$番目の地点 ($j=1,\ldots,n$) にどの程度重みを与えるかを表す量で，以下の形がよく用いられます．

w_{ij}=\exp\left(-\frac{\|s_i-s_j\|^2}{2h^2}\right)．

ただし，$h$はバンド幅と呼ばれる量で，どの程度遠い地点まで重みを与えるかを表すチューニングパラメータです．上の定義から，$s_i$と$s_j$が離れるにつれて$w_{ij}$が小さくなり，その減少度合いが$h$によって規定されていることがわかります．また$w_{ii}=1$ (同一地点の重みは1) であることもわかります．

以上から，(1)の推定方法の直感的なイメージとして，「$s_i$にある程度近い位置にあるデータは近似的に同じ回帰係数$\beta_i$をもつ回帰モデルに従う」ことを想定した局所的な最小二乗法を実行していると解釈することができます．

GWR推定量の形

(1)の推定量は解析的な解を与えることができます．
$y=(y_1,\ldots,y_n)$, $X=(x_1,\ldots,x_n)^t$, $W_i={\rm diag}(w_{i1},\ldots,w_{in})$と定義すると，

\widehat{\beta}_i=(X^tW_iX)^{-1}X^tW_iy  \tag{2}

となります．この量を各地点で計算することで空間的に異なる回帰係数を推定することができ，その違いを見ることで空間的な構造の違いを把握することができます．

バンド幅の選択

GWRを実行するためにはバンド幅$h$を適切に選択する必要があります．まずは，以下のように$h$の値によってバイアス(推定手法の正しさ)と分散(安定性)の間にトレードオフの関係があることに注意します．

小さい$h$: 局所的なデータのみを使うのでバイアスなく推定できるが，使えるデータが少ないため推定の分散が大きくなる．
大きい$h$: 広範囲のデータを使うため安定的に推定できる一方で，構造が異なるデータも含んでしまう可能性があるため推定にバイアスが生じる．

データから最適な$h$を決定する方法として以下のようなcross validation (CV)の基準を用います．

CV(h) = \sum_{i=1}^n\{y_i-x_i^t\widehat{\beta}_{i, -i}(h)\}^2．

ただし，$\widehat{\beta}_{i,-i}(h)$はバンド幅$h$のもとで$i$番目のデータ以外を用いた推定量を表します．すなわち，CV基準は**$y_i$以外のデータを用いて$y_i$を予測したときにどの程度予測が正確か**を表した量になります．いくつかの$h$の候補の中でCV基準を最小にするものを最適な$h$として採用します．

一般化

(1)の考え方は線形回帰モデルのみならず一般の統計モデルに拡張することができます．
今$f(y_i|x_i;\theta_i)$を一般の統計モデルとします（例えば一般化線形回帰モデルなど)．
ここで$\theta_i$は地点ごとに異なるパラメータで，(1)の$\beta_i$に相当するものです．
GWRと同じ考え方を用いて，$\theta_i$を以下のように推定することができます．

\widehat{\theta}_i={\rm argmax}_{\theta}\sum_{j=1}^n w_{ij}\log f(y_j|x_j;\theta)．

上記の目的関数は地点$i$周辺のデータのみで構成した局所的な対数尤度関数であると解釈することができます．さらに、$w_{ij}$に含まれるバンド幅$h$は以下のCV基準に基づき決定することができます．

CV(h) = -\sum_{i=1}^n \log f(y_i|x_i;\hat{\theta}_{i,-i}(h))．

ただし，$\hat{\theta}_{i,-i}(h)$はバンド幅$h$のもとで$i$番目のデータ以外を用いた推定量を表します．

Rでの実行

GWRはR上でspgwrパッケージを用いて簡単に実行することができます．基本的に使うのは以下の2つの関数です．

gwr.sel: CV規準により最適なバンド幅を選択する関数
gwr: バンド幅を与えたもとでGWRを実行する関数

詳細についてはspgwrパッケージのマニュアルをご覧ください．

試しにmlbenchパッケージに入っているBostonHousingデータに適当してみます．このデータは1970年代後半におけるボストンの国勢統計区ごとの住宅価格の中央値（被説明変数）と計13個の特徴量（説明変数）に関するデータで，機械学習や統計学におけるベンチマーク的なデータとして知られています．
データの詳細についてはこちらを参照してください．

library(spgwr)
library(mlbench)
data(BostonHousing2)

# データの用意
Y <- BostonHousing2$cmedv　   # 被説明変数
X <- subset(BostonHousing2, T, c(crim, zn, indus, chas, nox, rm, age, dis, rad, tax, ptratio, b, lstat))
X$chas <- as.numeric(X$chas) - 1
X <- as.matrix(X)   # 説明変数
Sp <- as.matrix( subset(BostonHousing2, T, c(lon, lat)) )      # 緯度・経度

# CVによる最適なバンド幅の選択
h <- gwr.sel(Y~X, coords=Sp)

# 最適なバンド幅を用いたGWR
fit <- gwr(Y~X, coords=Sp, bandwidth=h)

# 結果の図示
est <- fit$SDF$Xage
value <- (est - min(est)) / diff(range(est))
cs <- colorRamp( c("blue", "green", "yellow", "red"), space="rgb")
cols <- rgb( cs(value), maxColorValue=256 )
plot(Sp, col=cols, xlab="Latitude", ylab="Longitude", pch=20, cex=1)

age(1940年以前に建てられた住宅の割合)の回帰係数の推定値の空間分布は以下のようになります．
(青→黄→赤ほど値が高くなることを表します．)

この図から地点ごとに異なった回帰係数の値が推定されていることがわかります．

おわりに

GWRは簡単に実行ができて便利な方法ですが，以下のような弱点が知られています．

空間的に端の点や周囲のデータが少ない地点では数値的に不安定な結果になってしまう．(局所的に回帰モデルを当てはめるため．)
サンプル地点が数千以上の大規模空間データでは計算に非常に時間がかかってしまう．(各サンプル地点ごとに回帰モデルを推定する必要があるため．)

このような問題点を解決する方法論については近年研究が進んでおりますので，今後紹介していきたいと思います．また，空間異質性のもとでのより高度な回帰分析や理論については栗栖先生の記事もご参照ください．

株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospareまでお問い合わせください．