東京大学・株式会社Nospareの菅澤です．
今回は最近の研究成果「空間クラスタリングを用いた回帰分析」について紹介させていただきます．

はじめに: 地理的加重回帰(GWR)の弱点

空間データに対する標準的な回帰分析の方法として地理的加重回帰(GWR; Geographically Weighted Regression)があります．端的に言うと，GWRは（通常の線形回帰と異なり）空間的に回帰係数が変化する構造を考慮した手法です．詳しい解説については過去に書いたこちらの記事をご覧ください．

GWRは空間データを解析するための有効なツールですが，以下のような問題点が知られています．

数値的に不安定である．（局所的に回帰モデルを当てはめるため．)
サンプル地点が数万以上の大規模空間データでは計算コストが高い．（各サンプル地点ごとに回帰モデルを推定する必要があるため．）

１つ目の欠点は空間的に端の点など周囲のデータが少ない地点で顕著になることが知られています．近年では数万オーダー以上の大規模空間データが利用可能になっている一方で，サンプル地点は必ずしも空間的に均一に分布するとは限らない（例えば不動産のデータは都心部では密なのに対して郊外では疎な傾向にある）ため，上記の2つの問題点はデータ分析の場面でどちらも重要になってきます．

正則化を用いた新しい空間回帰の方法（Li and Sang, 2019)

提案手法を紹介する前に関連した既存研究を紹介します．
Li and Sang (2019)では，正則化の考え方を空間回帰分析に導入した方法を提案しています．GWRと同様に，以下のような回帰モデルを考えます．

y(s_i)=x(s_i)^T\beta(s_i)+\varepsilon(s_i), \ \ \ \  i=1,\ldots,n.

ここで$s_i$は位置情報（緯度・経度など）を表し，$y(s_i), x(s_i)$はそれぞれ被説明変数と説明変数ベクトルであり，$\beta(s_i)=(\beta_1(s_i),\ldots,\beta_p(s_i))^T$は空間的に変化する回帰係数です．また，$\varepsilon(s_i)$は空間的に独立な誤差項を表します．

GWRは局所的なウエイトを用いて$\beta(s_i)$を推定していましたが，Li and Sang (2019)は以下の目的関数を最小化することで$\beta(s_i)$を推定することを提案しました．

\frac1n\sum_{i=1}^n\{y(s_i)-x(s_i)^T\beta(s_i)\}^2 + \lambda\sum_{k=1}^p\sum_{i\sim j}|\beta_k(s_i)-\beta_k(s_j)|

ここで，$\lambda$はチューニングパラメータであり，$i\sim j$は2つの地点$s_i, s_j$が"隣接"していることを表します．データから隣接構造を決定する方法はいくつか考えられますが，基本的に「隣接している $=$ 空間的な距離が近い」というイメージです．

上の目的関数の第二項はL1正則化（特にfused lasso）と呼ばれる形になっており，$|\beta_k(s_i)-\beta_k(s_j)|$が部分的に正確に0になるような推定結果を返すデザインになっています．すなわち，上の目的関数を最小化することで，空間的に近い2つの係数$\beta_k(s_i), \beta_k(s_j)$は同じ値を取りやすい構造を推定することができます．これは，空間が有限個のクラスターに分割され，同じクラスター内では回帰係数の推定値が一定になることを意味しています．

論文中ではGWRに対する推定精度の改善が数値実験を通して示されています．一方で，この方法はサンプル数に応じて正則化項の数が増加していくため，大規模空間データに対する計算コストはそれなりに大きく，GWRの欠点2を解決できているとは言い難いです．また，上記の正則化による方法では暗に誤差項$\varepsilon(s_i)$の分散が空間的に一定であることを仮定しているため，データによってはその構造が制約的になり全体的な推定精度が低下してしまう恐れがあります．

提案手法: Spatially Clustered Regression (SCR)

今回開発した方法では，明示的な空間クラスタリング構造を統計モデルに導入することでGWRが持つ2つの欠点を同時に解決します．

以下ではより一般的な設定を考えます．$y_i$, $x_i$をそれぞれ被説明変数および説明変数ベクトルとし，$f(y_i|x_i;\theta_i)$を当てはめたい統計モデルとします．ここで$\theta_i$はモデルパラメータで空間的に変化することを許す設定になっています．例えばモデルとして$\phi(y_i;x_i^T\beta_i,\sigma^2)$を考えると前述のLi and Sang (2019)と同じ設定になります．

$\theta_i$を推定するために明示的なグルーピングパラメータ$g_i\in\{1,\ldots,G\}$($G$はグループ数）を導入します．
推定には以下の関数を最大化します．

Q(\theta, g)=\sum_{i=1}^n \log f(y_i|x_i;\theta_i) + \phi \sum_{i<j}\omega_{ij}I(g_i=g_j)

ただし，$\omega_{ij}$は地点$i$と地点$j$の重みを表します．（簡単なものとしては地点間の距離がある一定の値を下回る場合は$\omega_{ij}=1$，そうでない場合は$\omega_{ij}=0$とする方法が考えられます．）
また，$\phi$はグルーピングの空間的な関連の強さを表すチューニングパラメータで，以下では簡単のため$\phi=1$としておきます．

第二項は「空間的に近い地点同士が同じグループに含まれることで目的関数の値を増加させる」役割を果たしていることがわかります．すなわち$Q(\theta, g)$を最大化することで，空間的に近い地点同時が同じグループに属するようなクラスタリング構造を得ることができます．また，この関数の最大化はk-meansクラスタリングの繰り返しアルゴリズムと似た方法で簡単に実行することができます．

ちなみに，第二項で登場する（ある種の）罰則項はSugasawa (2020)で提案されたもので，離散値をとる空間モデルとして知られているPottsモデルに由来しています．

提案手法を実際に使うためにはクラスター数（$G$）を適切に決定する必要がありますが，今回は単純にBICを用いて選択する方法を採用します．

数値実験による比較

以下のデータ生成過程を考えます．

y_i=\beta_{0i}+\beta_{1i} x_{i1}+\beta_{2i}x_{2i}+\varepsilon_i, \ \ \ \varepsilon_i\sim N(0, \sigma_i^2), \ \ \ i=1,\ldots,1000.

回帰係数($\beta_{0i}, \beta_{1i}, \beta_{2i}$)および誤差分散($\sigma_i^2$)の真の構造が

（１）空間的なクラスタリング構造がある
（２）空間的に滑らかに変化する

の２つのケースを考えて提案手法と既存手法のパフォーマンスを比較しました．(2)のケースでは真の構造がクラスタリング構造を有していないため，提案手法で仮定されている構造は真の構造をカバーできていないのですが，近似的なモデルとしてどのくらい適合的なパフォーマンスを示すのかを調べる目的のシナリオです．

以下の方法を用いて回帰係数の推定を行います．

GWR: 地理的加重回帰 (バンド幅はcross validationで選択する)
SHP: Li and Sang (2019)の方法 (正則化パラメータはBICで選択する)
SCR: 提案手法 (クラスター数はBICで選択する)

結果

下の図は$\beta_{0i}$ (左)，$\beta_{1i}$ (中央)，$\beta_{2i}$ (右)の真の値（True）とそれぞれの手法の推定結果を表しています．

（１）空間的なクラスタリング構造があるケースの結果

クラスタリング構造をもつSHPとSCRは真の構造をうまく推定できていることが観察できます．一方でGWRは所々推定値の変化が激しい箇所が出てきています．

（２）空間的に滑らかに変化するケースの結果

真の構造がクラスタリング構造を有していなくても，SCRは適切に真の構造を推定できていることがわかります．カラクリとしては，このケースにおいてBICが選択したクラスター数は30であり，回帰係数が空間的に滑らかに変化する場合においても大きいクラスター数を用いることで適合的に推定を実行できていることがわかります．

計算時間の比較

今度は提案手法の計算時間的な効率性について検証してみます．サンプル数を1000, 3000, 5000, 10000, 20000と増やしていき，計算時間がどのよう推移していくかを確認しました．（SHPは1000個のサンプルのケースにおいて既に他と比べてそれなりの計算時間がかかっていたため，今回の比較からは除外しています．）また，以下の結果はチューニングパラメータの選択も込みの計算時間です．

結果は以下のようになりました．
(SFCRについては本記事では紹介しておりませんが，SCRの亜種となる手法です．)

この結果を見ると，サンプル数が20000の場合，GWRは16分程度かかるのに対しSCRは1分程度しかかからないことがわかります．

Rコード

提案手法を実装したコードや上記の数値実験で用いた擬似データがGitHubからダウンロードできます．

提案手法はSCR-function.Rというファイルに実装されています．以下は上記の数値実験の結果を再現するためのコードです．

# 必要なパッケージのインストール
library(spgwr)   
library(glmnet) 
library(ape)   

# 提案手法のインストール
source("SCR-function.R")

# 擬似データのロード
scenario <- 1    # 1(シナリオ1)か2(シナリオ2)を選ぶ
load(paste0("simdata", scenario, ".RData"))

# GWRの当てはめ
b.opt <- gwr.sel(Y~X, coords=Sp, verbose=F)
fit <- gwr(Y~X, coords=Sp, bandwidth=b.opt)
est.GWR <- cbind(fit$SDF$`(Intercept)`, fit$SDF$`Xx1`, fit$SDF$`Xx2`)

## SCRの当てはめ
# 最適なクラスター数の選択
G.set <- seq(5, 30, by=5)
IC <- SCR.select(Y, X, W, Sp, G.set=G.set, Phi=1, print=F, family="gaussian")
hG <- IC$select[2]     

# 最適なクラスター数のもとでの推定 
fit.SCR <- SCR(Y, X, W, Sp, G=hG, Phi=1)   
est.SCR <- fit.SCR$sBeta


## SHPの当てはめ
# minimum spanning treeを用いた正則化項の構成
n <- length(Y)
MST <- mst(dist(Sp))    
H <- c()
for(i in 1:n){
  for(j in 1:n){
    if(i<j){
      if(MST[i,j]==1){
        h <- rep(0, n)
        h[i] <- 1
        h[j] <- -1
        H <- rbind(H, h)
      }
    }
  }
}
HH <- rbind(H, rep(1/n, n))
invH <- solve(HH)
XX <- cbind(diag(n), diag(X[,1]), diag(X[,2]))
invHH <- matrix(0, 3*n, 3*n)
for(k in 1:3){
  sub <- (1+n*(k-1)):(n*k)
  invHH[sub, sub] <- invH
}
XH <- XX%*%invHH
pen <- rep(c(rep(1, n-1), 0), 3)

# 複数の正則化パラメータのもとでの当てはめ
fit.SHP <- glmnet(x=XH, y=Y, family="gaussian", intercept=F, penalty.factor=pen)

# 最適な正則化パラメータの選択
BIC <- c()
for(j in 1:100){
  mu <- as.vector(XH%*%as.vector(fit.SHP$beta[,j]))
  sig <- sqrt(mean((Y-mu)^2))
  BIC[j] <- -2*sum(dnorm(Y, mu, sig, log=T)) + log(n)*fit.SHP$df[j]
}
opt <- which.min(BIC)
est <- as.vector(invHH%*%as.vector(fit.SHP$beta[,opt]))
est.SHP <- matrix(est, n, 3)


## 結果の可視化
# 空間プロット用の関数
SPlot <- function(value, ran, xlim=NULL, title=""){
  value <- (value - ran[1]) / diff(ran)
  cs <- colorRamp( c("blue", "green", "yellow", "red"), space="rgb")
  cols <- rgb( cs(value), maxColorValue=256 )
  ran[1] <- floor(ran[1])
  ran[2] <- ceiling(ran[2])
  plot(Sp, col=cols, xlab="Latitude", ylab="Longitude", xlim=c(-1, 1.5), main=title)
  cs <- colorRamp( c("blue", "green", "yellow", "red"), space="rgb")
  cols <- rgb(cs(0:1000/1000), maxColorValue=256)
  rect(1.1, seq(0,2,length=1001), 1.3, seq(0,2,length=1001), col=cols, border=cols)
  tx <- seq(ran[1], ran[2], length=5)
  text(x=1.4, y=seq(0, 2, length=5), tx, cex=0.7)
  yy <- seq(0, 2, length=5)
  for (i in 1:5){
    segments(1.1, yy[i], 1.3, yy[i], col="white")
  }
}

# プロット
par(mfcol=c(4,3))
for(k in 1:3){
  ran <- range(Beta[,k], est.GWR[,k], est.SHP[,k], est.SCR[,k])
    ran[1] <- floor(ran[1])
    ran[2] <- ceiling(ran[2])
    SPlot(Beta[,k], ran, title="True")
    SPlot(est.GWR[,k], ran, title="GWR")
    SPlot(est.SHP[,k], ran, title="SHP")
    SPlot(est.SCR[,k], ran, title="SCR")
}

まとめ

本研究では空間データの回帰分析のための新しい方法論について紹介しました．GWRに替わる効果的な手法として様々なデータ（不動産データ，犯罪データ，疫学データなど）で使える方法ではないかと思います．

詳しい内容については以下の論文をご参照ください．

Sugasawa, S. and Murakami, D. (2020). Spatially clustered regression. arXiv:2011.01493. (投稿中)

今回の内容に関連するお問い合わせにつきましては，お気軽に菅澤までご連絡ください．

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