空間的異質性を考慮したノンパラメトリック回帰モデル

東京工業大学/株式会社Nospare リサーチャーの栗栖です．
この記事では最近の研究成果「空間的異質性を考慮したノンパラメトリック回帰に関する論文」(Kurisu(2021)) について紹介します．

モチベーションと 空間的異質性

分析したいデータが空間内のある領域で観測され，異なる地点のデータが互いに何らかの関係性をもつと考えられる場合，そのようなデータは空間データと呼ばれます．例えばある地点における気温や風速，降水量などの気候データ，$\text{PM}_{2.5}$や一酸化炭素といった化学物質の濃度など環境データ，ある地域における土壌調査や植生調査で得られたデータ，地価や住宅価格のデータなどは空間データとして扱うことができます．また空間データを扱う統計学の分野は空間統計と呼ばれます．空間統計では，しばしばデータの観測領域においてその性質が変化しないための仮定として，データ分析の際に用いる空間モデルに対して定常性と呼ばれる性質が仮定されることがあります．一方，空間データの観測範囲に異なる気候の地域が複数含まれていたり，異なる経済・人口規模の地域が含まれていたりすると，観測領域内においてデータが異なる空間的性質をもつ場合があります．このような場合データは空間的異質性 (または空間的非定常性)をもつといい，定常性の下でのデータ分析手法とは異なる分析手法が必要になります．

この研究では空間的異質性をもつデータに対する分析手法として，空間的な非定常性をもつノンパラメトリックな回帰モデルを新たに提案しました．特に，このモデルの (ノンパラメトリックな) 回帰関数を推定するために，伝統的な統計手法であるノンパラメトリックカーネル回帰推定量がどのような性質を持つかを理論的に調べました．

モデル

以下ではこの研究で提案したモデルについて紹介したいと思います．

モデルで扱う変数の設定

$Y_{A_n}(\boldsymbol{s}) \in \mathbb{R}$ : 地点 $\boldsymbol{s}$ での被説明変数，
$\boldsymbol{X}_{A_n}(\boldsymbol{s}) = \left(X^{(1)}_{A_n}(\boldsymbol{s}), \cdots, X^{(p)}_{A_n}(\boldsymbol{s})\right)' \in \mathbb{R}^{p}$ : 地点 $\boldsymbol{s}$ での $p$ 変量の説明変数．
添え字 $A_n$ は $Y_{A_n}(\boldsymbol{s})$, $\boldsymbol{X}_{A_n}(\boldsymbol{s})$ がともに観測領域内で空間的異質性をもつことを明示するために付けています．

データ観測の枠組みの設定

空間内の領域 $R_{n} = [0,A_n]^d \subset \mathbb{R}^d$ 内の $n$ 地点 $\boldsymbol{s}_1$, $\cdots$, $\boldsymbol{s}_n$ でデータを観測 (空間dデータでは特に $d=2$)．
すなわち，$(Y_{A_n}(\boldsymbol{s}_{1}), \boldsymbol{X}'_{A_n}(\boldsymbol{s}_{1}))'$, $\cdots$, $(Y_{A_n}(\boldsymbol{s}_{n}), \boldsymbol{X}'_{A_n}(\boldsymbol{s}_{n}))'$ が観測値(データ)です¹．
これらを使って以下で紹介するモデルの回帰関数 $m$ を推定することがこの研究の目標です．
ただし，観測地点 $\boldsymbol{s}_j$ は観測領域 $R_n$ 内で不等間隔に分布している場合を想定します．

データ観測の枠組みについて，観測地点は規則的に配置された格子点上に分布する，という流儀もあります．ただし，各格子点に画素が与えられた画像データのような場合を除き，空間データは物理的制約により観測地点が不等間隔になる場合が多いため，この研究では格子点上でデータが観測される場合も特別な場合として扱える上記の設定で考えることにします．

[例] 下の図はアメリカ全土にある気象台の配置図で，それぞれの点が気象台の位置に対応しています²．この図を見ると気象台は比較的密に配置されている地域もあればそうでない地域もあり，不等間隔に配置されていることが見て取れます．またアメリカは国土が広く，乾燥した地域や雨・雪の多い地域，寒冷地や温暖な地域が含まれるため，州・地域ごとの降水量などの気象データは空間的異質性を持つと考えられます．

一般モデル

この論文ではまず $Y_{A_n}(\boldsymbol{s})$ と $\boldsymbol{X}_{A_n}(\boldsymbol{s})$の間に次のような関係がある場合を考えます:

\begin{align}
Y_{A_n}(\boldsymbol{s}) = m\left({\boldsymbol{s} \over A_n}, \boldsymbol{X}_{A_n}(\boldsymbol{s})\right) + \varepsilon_{A_n}(\boldsymbol{s}).
\end{align}

ここで，$\varepsilon_{A_n}(\boldsymbol{s})$は観測誤差を表現していて，$\boldsymbol{X}_{A_n}(\boldsymbol{s})$ で条件付けた時の期待値は $0$，即ち $E[\varepsilon_{A_n}(\boldsymbol{s})|\boldsymbol{X}_{A_n}(\boldsymbol{s})] = 0$と仮定します．
このとき，関数 $m(\boldsymbol{u}, \boldsymbol{x})$ を (ノンパラメトリックな) 回帰関数と呼び ($\boldsymbol{u} \in [0,1]^{d}$, $\boldsymbol{x} = (x_{1},\cdots,x_{p})' \in \mathbb{R}^{p}$)，データから $m$ を推定することが目標です．特に線形回帰モデル (パラメトリックモデル) の場合は

Y_{A_n}(\boldsymbol{s}) =  \beta_{0} + \beta_{1}X^{(1)}_{A_n}(\boldsymbol{s}) + \cdots + \beta_{p}X^{(p)}_{A_n}(\boldsymbol{s}) + \varepsilon_{A_n}(\boldsymbol{s})

と表現することができ，この場合，$m(\boldsymbol{u}, \boldsymbol{x}) = m(\boldsymbol{x}) = \beta_{0} + \beta_{1}x_{1} + \cdots + \beta_{p}x_{p}$ であり，$m$ は空間内の地点 $\boldsymbol{s}$ には依存しません．一方で一般モデルでは回帰関数が($\boldsymbol{u}$にも依存するので)各地点ごとに変化する場合も考慮できる点が重要です．論文中では以下で定義されるカーネル推定量 $\hat{m}(\boldsymbol{u},\boldsymbol{x})$ を用いて $m(\boldsymbol{u},\boldsymbol{x})$ を推定する方法を提案し，その理論的性質を調べています．

\hat{m}(\boldsymbol{u},\boldsymbol{x}) = {\sum_{j=1}^{n}K_{h}\left(\boldsymbol{u} - {\boldsymbol{s}_{j} \over A_{n}},\boldsymbol{x} - \boldsymbol{X}_{A_{n}}(\boldsymbol{s}_{j})\right)Y_{A_{n}}(\boldsymbol{s}_{j}) \over \sum_{j=1}^{n}K_{h}\left(\boldsymbol{u} - {\boldsymbol{s}_{j} \over A_{n}},\boldsymbol{x} - \boldsymbol{X}_{A_{n}}(\boldsymbol{s}_{j})\right)},

ただし，$K_{h}(\boldsymbol{u},\boldsymbol{x}) = K(\boldsymbol{u}/h,\boldsymbol{x}/h)$，
$h$ はバンド幅と呼ばれ，$n \to \infty$ で $h \to 0$ となる定数列，
$K:\mathbb{R}^{d+p} \to \mathbb{R}$ はカーネル関数と呼ばれ，この論文では原点対称でコンパクトなサポートを持ち，リプシッツ連続な関数としています.

特に $n \to \infty$ のとき，$h$ を適切にとることで $\hat{m}(\boldsymbol{u},\boldsymbol{x})$ が $m(\boldsymbol{u},\boldsymbol{x})$ に (確率)収束することが示すことができます．

加法モデル

一般モデルは柔軟なモデリングが可能である一方，説明変数の次元 $p$ が大きくなるにつれて回帰関数 $m(\boldsymbol{u}, \boldsymbol{x})$ の推定量 $\hat{m}(\boldsymbol{u}, \boldsymbol{x})$ が $m(\boldsymbol{u}, \boldsymbol{x})$ に収束するスピードが遅くなってしまうという問題 (次元の呪い)があります．この問題に対処するため，なるべく柔軟性を損なわず，なおかつ次元の呪いの影響を受けないモデルとして，この論文では次のような加法モデルも考えています³．

Y_{A_n}(\boldsymbol{s}) =m_{0}\left({\boldsymbol{s} \over A_n}\right) + m_{1}\left({\boldsymbol{s} \over A_n}, X^{(1)}_{A_n}(\boldsymbol{s})\right) + \cdots + m_{p}\left({\boldsymbol{s} \over A_n}, X^{(p)}_{A_n}(\boldsymbol{s})\right) + \varepsilon_{A_n}(\boldsymbol{s}).

このモデルは一般モデルで $m(\boldsymbol{u}, \boldsymbol{x}) = m_{0}(\boldsymbol{u}) + m_{1}(\boldsymbol{u}, x_{1}) + \cdots + m_{p}(\boldsymbol{u}, x_{p})$ となる場合に対応しています．回帰関数 $m(\boldsymbol{u},\boldsymbol{x})$ が説明変数の各成分ごとにモデル化された関数 $m_{j}(\boldsymbol{u}, x_{j})$ (または $m_{0}(\boldsymbol{u})$) の和の形で書かれていて，各 $m_{j}(\boldsymbol{u},x_{j})$ は依然としてノンパラメトリックな関数であることがポイントです．加法モデルは以下のような (セミパラメトリック) 線形回帰モデルも含みます．

Y_{A_n}(\boldsymbol{s}) = \beta_{0}\left({\boldsymbol{s} \over A_n}\right) + \beta_{1}\left({\boldsymbol{s} \over A_n}\right)X^{(1)}_{A_n}(\boldsymbol{s}) + \cdots + \beta_{p}\left({\boldsymbol{s} \over A_n}\right)X^{(p)}_{A_n}(\boldsymbol{s}) + \varepsilon_{A_n}(\boldsymbol{s}).

上記のモデルは特に空間係数変化モデル (spatially varying coefficient model) と呼ばれ，通常の線形回帰モデルの係数が各地点ごとに変化するモデルになっており，係数が一定の線形回帰モデルにくらべてより柔軟なモデリングが可能になっています．論文中では一般モデルと同様にカーネル推定量を用いて $\hat{m}_{0}(\boldsymbol{u})$, $\hat{m}_{1}(\boldsymbol{u},\boldsymbol{x})$,$\cdots$, $\hat{m}_{p}(\boldsymbol{u},\boldsymbol{x})$ を構成する方法を提案し，推定量の性質を調べています．

応用例

この研究で提案した方法の具体的な応用例としては，
[例1] 気象データ(降水量などの気象データを用いた地域特有の性質(空間的異質性)の検出，Matsuda and Yajima(2018))
[例2] 環境データ(ある沿岸水域における化学物質や有機物の濃度が水質に与える影響分析，Hallin et al.(2004))
[例3] 経済データ(地価や住宅価格などの不動産データ分析，Gelfand et al.(2003))
などが挙げられます．

モデルの特徴

ここまでに紹介したモデルの特徴についてまとめておきます．

(1) データの観測地点が不等間隔でもOK．
$\to$ 空間データは物理的制約により観測地点が不等間隔になる場合が多い．
(例) 降水量・気温・風速・大気の状態などの観測所，土壌調査のサンプル採取地点，
ある地域における植生の調査地点，地価データ，住宅価格データ etc.

(2) 回帰関数が空間内の各地点に依存して異なる場合(空間的異質性がある場合)にも対応できる．
$\to$ 空間データは各地点あるいは各地域ごとに変化する特徴を持っていると考えたほうが自然な場合が多い．
(例) 標高が高い or 低い, 都心 or 郊外 etc.

(3) 被説明変数と説明変数の関係の間に非線形な関係がある場合にも対応できる．
$\to$ 線形回帰モデルよりも柔軟なモデリングが可能．

(4) 加法モデルを考えることでモデルの柔軟性を保ちつつ，次元の呪いの影響を取り除くことができる．

理論的内容について(おまけ)

モデルの導入の際には空間的異質性をもつデータの具体的な定式化には触れませんでしたが，論文の中では局所定常性(local stationarity) と呼ばれる考え方を用いて空間的な非定常性を表現できるようにしています．局所定常性とは直感的には，データ自身は非定常だが，ある微小な領域ではそのデータは定常な空間モデルで近似できる，という性質です．この性質により，定常性の下での理論解析が空間非定常な場合にも拡張できるようになりました．

さらに論文では一般モデルに対して，回帰関数 $m$ の カーネル推定量 $\hat{m}$ の (あるコンパクト集合上での一様な) 一致性と各点での漸近正規性を示しています．また加法モデルに対しても，カーネル推定量と smooth backfitting method (Mammen et al. (1999)) と呼ばれる方法を組み合わせて各 $m_{j}$ が推定可能であることを示し，一般モデルの場合と同様に各 $m_j$ のカーネル推定量 $\hat{m}_j$ の (あるコンパクト集合上での一様な) 一致性と各点での漸近正規性を示しています．

局所定常性の数学的な定式化やノンパラメトリックカーネル回帰推定量の構成とその理論的な解析の詳細については Kurisu(2021)，その他の記事の中で紹介した論文については以下の参考文献をご参照ください．

まとめ

この記事では空間的な性質が地域によって異なる場合におけるノンパラメトリック回帰分析について最近の研究成果を紹介しました．株式会社Nospareでは空間統計に限らず，統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータ分析につきましては株式会社Nospareまでお問い合わせください．

参考文献
[1] Gelfand, A.E., Kim, H.-J., Sirmans, C.F. and Banerjee, S. (2003) Spatial modeling with spatially varying coefficient processes. Journal of the American Statistical Association 98, 387-396.
[2] Hallin, M., Lu, Z. and Yu, K. (2009). Local linear spatial quantile regression. Bernoulli 15, 659-686.
[3] Kurisu, D. (2022) Nonparametric regression for locally stationary random fields under stochastic sampling design. arXiv:2005.0637 Bernoulli 28, 1250-1275.
[4] Mammen, E., Linton, O. and Nielsen, J. (1999) The existence and asymptotic properties of a backfitting projection algorithm under weak conditions. Annals of Statistics 27, 1443-1490.
[5] Matsuda, Y. and Yajima, Y. (2018) Locally stationary spatio-temporal processes. Japanese Journal of Statistics and Data Science 1, 41-57

1. 記事では $R_{n}$ は長方形としていますが，多角形や楕円形のような凸集合，星形のような非凸集合にも $R_{n}$ の定義を拡張することができます. ↩

2. 出典：Institute of Mathematics Applied to Geosciences (IMAGe) https://www.image.ucar.edu/Data/US.monthly.met/ ↩

3. 一般モデルでは $\hat{m}(\boldsymbol{u},\boldsymbol{x})$ の収束レートが $O_{p}\left(\sqrt{\log n \over nh^{d+p}}\right)$ であるのに対し，加法モデルでは， $m_{0}$ については $O_{p}\left(\sqrt{\log n \over nh^{d}}\right)$, $m_{1},\cdots,m_{p}$ にたいしてはそれぞれ $O_{p}\left(\sqrt{\log n \over nh^{d+1}}\right)$ であることがわかります．特に $p \geq 2$ の場合は加法モデルの方が収束レートが早くなることがわかります． ↩