局所多項式回帰による空間データ分析（１）

はじめに

東京大学/株式会社Nospare リサーチャーの栗栖です．
この記事では最近の研究成果「空間データに対する局所多項式トレンド回帰」(Kurisu and Matsuda(2023+))について紹介します．本研究の目的は，ある地域内の様々な地点$\boldsymbol{X}_1,\dots,\boldsymbol{X}_n$で観測されるデータ$Y_1=Y(\boldsymbol{X}_1),\dots,Y_n = Y(\boldsymbol{X}_n)$を用いてその地域内の空間データのトレンド(空間トレンド)を推定することです．

[例] 下の図はアメリカ全土にある気象台の配置図で，それぞれの点が気象台の位置に対応しています¹．各気象台では毎月の降水量が記録されています．例えばある月のアメリカ国内の様々な地域において降水量がどのように変化しているかに興味があるとします．このとき，地点$\boldsymbol{X}_i$において実際に記録された降水量$Y_i=Y(\boldsymbol{X}_i)$は

Y_i = m(\boldsymbol{X}_i) + v_i

という構造をしていると仮定します．ここで$m(x)$はアメリカ国内の様々な地域における降水量の変化を説明する関数で，$v_i$は各地点における観測誤差(観測地点ごとの気象条件による降水量の観測誤差や計測誤差等を含めたもの)です．この例では，アメリカ国内の各気象台とそこで観測された降水量のデータ$(Y_i,\boldsymbol{X}_i)_{i=1}^n$ (ここではサンプルサイズ$n$は気象台の数)を用いて関数$m$(空間トレンド)を推定することが目的となります．

モデル

2次元の空間 $R = [-A,A]^2$ 内の地点 $\boldsymbol{X}_i=(X_{1,i},X_{2,i})'$, $i=1,\dots,n$ において空間データ $Y_i = Y(\boldsymbol{X}_i)$, $i=1,\dots,n$ を観測する状況を考えます．このとき，データ $(Y_i, \boldsymbol{X}_i)$, $i = 1,\dots,n$ が以下のノンパラメトリックな回帰モデルに従うとします：

Y_i = m\left({\boldsymbol{X}_i \over A}\right) + e(\boldsymbol{X}_i) + \varepsilon_i,\ i=1,\dots,n

ここで，$\boldsymbol{e}=\{e(\boldsymbol{x}): \boldsymbol{x} \in \mathbb{R}^2\}$ は空間相関をもつ誤差項で，データ$Y_i$ の空間データとしての構造を表現する部分です．特に $E[e(\boldsymbol{X}_i)|\boldsymbol{X}_i=x]=0$ を満たすとします．また $\varepsilon_1,\dots, \varepsilon_n$ は各観測地点における独立同分布な観測誤差で，$E[\varepsilon_i]=0$ を満たすとします．空間相関をもつ誤差項がなければ上記のモデルは通常の独立同分布なデータに対するノンパラメトリック回帰モデルとなります．この意味で，上記のモデルは独立なデータに対するノンパラメトリック回帰分析の枠組みを自然に空間データに拡張したモデルとなっています．

以下ではデータから回帰関数 $m(x) = E[Y_i|\boldsymbol{X}_i/A = x]$, $x = (x_1,x_2) \in [-1,1]^2$ を推定するという問題を考えましょう．推定にあたり，上記のモデルに対して以下を仮定します．

(A1:SLP)
[1] $\boldsymbol{X}_i$は $R$上で定義された密度関数 $f$ をもち，$f$ は $x$ のある近傍 $N_x$ で連続かつ $f(x)>0$．
[2] $\{\boldsymbol{X}_i\}_{i=1}^{n}$，$\boldsymbol{e}=\{e(x): x \in \mathbb{R}^2\}$，$\{\varepsilon_i\}_{i=1}^{n}$ は独立．　
[3] $E[e^2(x)]=E[\varepsilon_i^2]=1$．ある整数 $q \geq 5$ に対して $E[|e(x)|^q], E[|\varepsilon_1|^q]<\infty$.
[4] $\sigma_e(x)=E[e(0)e(x)]$ として，$\int_{\mathbb{R}^2} |\sigma_e(x)|d x<\infty$．
[5] $m$ は $N_x$ で $p+1$ 回連続偏微分可能．

上記の条件[4]における $\sigma_e(\boldsymbol{x})$ は空間相関をもつ誤差項の相関関数に対応しています．
また $x \in \mathbb{R}^2$ で $p+1$ 回連続偏微分可能な関数 $g: \mathbb{R}^2 \to \mathbb{R}$ に対して，

\partial_{j_1 \dots j_L} g(x) = {\partial  \over \partial x_{j_1} \dots \partial x_{j_L}}g(x),\ 1\leq j_1 \leq \dots \leq j_L \leq 2,\ L=1,\dots, p+1

とします．

局所多項式回帰

回帰関数 $m$ の推定にあたり，以下の最適化問題を考えましょう．

\begin{align*}
\hat{\beta}(x) &= (\hat{\beta}_0, \hat{\beta}_{1},\hat{\beta}_{2},\hat{\beta}_{11},\hat{\beta}_{12},\hat{\beta}_{22}, \hat{\beta}_{111},\dots, \hat{\beta}_{222},\dots, \hat{\beta}_{1\dots 1}, \dots, \hat{\beta}_{2 \dots 2})\\ 
&= \underset{\beta \in \mathbb{R}^{D}}{\text{argmin}}\!\sum_{i=1}^n \!\!\left(\! Y_i \!-\! \beta_0 \!-\! \sum_{L=1}^{p}\sum_{1 \leq j_1 \leq \dots \leq j_L \leq 2}\!\!\!\!\!\!\!\!\!\!\! \beta_{j_1 \dots j_L}\! \prod_{\ell=1}^{L}\left({X_{j_\ell, i} \over A} - x_{j_\ell}\right)\!\!\right)^2 \!\!\! K \!\! \left(\!{\boldsymbol{X}_i/A - x \over h} \!\right),\\
K\left({\boldsymbol{X}_i/A - x \over h}\right) &= K\left({X_{1,i}/A - x_1 \over h},{X_{2,i}/A - x_2 \over h}\right).
\end{align*}

ここで，

D = 1 + \#\{(j_1,\dots, j_L): 1\leq j_1 \leq \dots \leq j_L\leq 2, 1\leq L\leq p\},

$h = h_{n} \to 0$, $n \to \infty$ はバンド幅，$K:\mathbb{R}^2 \to \mathbb{R}$ はカーネル関数と呼ばれ，以下を満たすとします．

(A2:SLP)
[1] $\int K(u)du=1$
[2] $K(u)$ は有界かつ $[-1,1]^2$ に台をもつ ($K(u) = 0$, $\|u\| >1$)
[3] $\tilde{u} = (u_{j_1\dots j_L})_{1\leq j_1 \leq \dots \leq j_L\leq 2, 1 \leq L \leq p}$, $\check{u} = (1, \tilde{u})'$ として，$S=\int \check{u}\check{u}'K(u)du$ は正則(逆行列をもつ)

上記の最適化問題は $\boldsymbol{X}_i/A = x \in [-1,1]^2$ の近傍で「局所」的に $Y_i$ を「多項式関数」に「回帰」していると解釈できるので，「局所多項式回帰」と呼ばれます．また得られる解は「局所多項式推定量」(local polynomial estimator, LP 推定量) と呼ばれることもあります．特にこの方法では，$\boldsymbol{X}_i/A = x$ の近傍で

\begin{align*}
Y_i &\approx m(\boldsymbol{X}_i/A) \approx m(x) + \sum_{L=1}^{p}\sum_{1 \leq j_1 \leq \dots \leq j_L\leq 2}{\partial_{j_1 \dots j_L} m(x) \over s_{j_1\dots j_L}!}\prod_{\ell=1}^{L}\left({X_{j_\ell, i} \over A} - x_{j_\ell}\right)
\end{align*}

と近似できることを利用して $m(x)$ (とその偏微分 $\partial_{j_1\dots j_L} m(x)$) を推定しています．ここで，

\begin{align*}
s_{j_1\dots j_L}! &= \prod_{k=1}^{2}s_{j_1\dots j_Lk}!,\\ 
s_{j_1\dots j_Lk} &= \#\{j_\ell: j_\ell = k, 1\leq \ell \leq L\}.
\end{align*}

上記の最適化問題を解くと，$m(x)$ の LP推定量が以下で与えられることが分かります．

\begin{align*}
\hat{\beta}(x) &= \left(\hat{m}(x),\partial_1 \hat{m}(x) ,\partial_2 \hat{m}(x), {\partial_{11}\hat{m}(x) \over s_{11}!},{\partial_{12}\hat{m}(x) \over s_{12}!},{\partial_{22}\hat{m}(x) \over s_{22}!} \right.,\\ 
&\left. \quad {\quad \partial_{111}\hat{m}(x) \over s_{111}!},\dots,{\partial_{222}\hat{m}(x) \over s_{222}!},\dots,{\partial_{1\dots 1}\hat{m}(x) \over s_{1\dots1}!},\dots,{\partial_{2 \dots 2}\hat{m}(x) \over s_{2 \dots 2}!}\right) \\ 
&=　\left(\sum_{i=1}^{n}K\left({\boldsymbol{X}_i/A - x \over h}\right)\tilde{\boldsymbol{X}}_i\tilde{\boldsymbol{X}}'_i \right)^{-1}\sum_{i=1}^{n}K\left({\boldsymbol{X}_i/A - x \over h}\right)\tilde{\boldsymbol{X}}_iY_i 
\end{align*}

ここで $\tilde{\boldsymbol{X}}_i = (1, \check{\boldsymbol{X}}_i)'$, $\check{\boldsymbol{X}}'_i = ((\boldsymbol{X}_i/A - x)'_1,\dots, (\boldsymbol{X}_i/A - x)'_p)'$,

(\boldsymbol{X}_i/A - x)_L = \left(\prod_{\ell=1}^{L} \left({X_{j_\ell,i} \over A} - x_{j_\ell}\right)\right)_{1 \leq j_1 \leq \dots \leq j_L \leq 2}

です．このとき，$\boldsymbol{X}_i/A = x$ の近傍で　

\begin{align*}
Y_i &\approx m(\boldsymbol{X}_i/A) \approx \tilde{\boldsymbol{X}}'_iM(x),\\
M(x) &= \left(m(x), \check{M}(x)\right)',\ \check{M}(x) = \left({\partial_{j_1 \dots j_L} m(x) \over s_{j_1\dots j_L}!}\right)'_{1 \leq j_1 \leq \dots \leq j_L\leq 2, 1 \leq L \leq p}
\end{align*}

と書けることに注意しておきます．従って，$\hat{\beta}$ は $M(x)$ の推定量とみることができます．

LP推定量の漸近正規性

最後にLP推定量の漸近的性質について紹介しておきます．まず(A1:SLP), (A2:SLP) に加えて以下を仮定します．

(A3:SLP)
[1] $n \to \infty$ として，$h \to 0$, $nh^2 \to \infty$．
[2] $n \to \infty$ として，$A/n \to \kappa \in [0,\infty)$，$Ah^{2+2p} \to \infty$，$Ah^{4+2p} \to c \in [0,\infty)$．

$H = \text{diag}(1,h,h, h^2,h^2,h^2, \dots, h^p, \dots, h^p) \in \mathbb{R}^D$ とします．
仮定(A1:SLP)-(A3:SLP)と適当な条件の下で，$n \to \infty$ として以下が成り立ちます．

\begin{align*}
&\sqrt{Ah^2}\left(H\left(\hat{\beta}(x) - M(x)\right)- S^{-1}B_p M_{p,n}(x)\right)\\ 
&\stackrel{d}{\to}N
\left( \boldsymbol{0}, \left\{{2\kappa \over f(x)} + \int_{\mathbb{R}^2}\sigma_e(x)dx\right\}S^{-1} V S^{-1}
\right).
\end{align*}

ここで，$V = \int \check{u}\check{u}'K^2(u)du$,

\begin{align*}
B_p &= \int \check{u}(u)'_{p+1}K(u)du,\ (u)_{p+1} = \left(\prod_{\ell=1}^{p+1}u_{j_\ell}\right)'_{1 \leq j_1 \leq \dots \leq j_{p+1} \leq 2},\\
M_{p,n}(x) &= \left({\partial_{j_1 \dots j_{p+1}} m(x) \over s_{j_1\dots j_{p+1}}!}h^{p+1}\right)'_{1 \leq j_1 \leq \dots \leq j_{p+1}\leq 2}.
\end{align*}

上記の結果において，推定量の漸近分布として現れる正規分布の分散の$\int_{\mathbb{R}^2}\sigma_e(x)dx$ の部分を形式的に $0$ とすると，独立同分布のデータに対するLP推定量の漸近的な分布に一致します．またデータが空間相関をもつ場合はLP推定量の収束レートが $\sqrt{Ah^2}$ となる点も特徴です (独立データの場合のLP推定量の収束レートは $\sqrt{nh^2}$)．これは，空間データ $Y_i$ の空間相関をとらえるには観測領域 $R = [-A,A]^2$ が拡大する必要があり，$R$ 内のデータ数 $n$ が増えても空間相関に関してデータがもつ情報が増えないためです．

次回の記事では，回帰関数の信頼区間の構成にあたり必要となる推定量の漸近分散の推定方法について紹介します．またこの記事では推定量の各点での漸近正規性について解説しましたが，次回の記事ではある集合上での推定量の一様収束レートについても解説する予定です．

まとめ

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

参考文献
[1] Kurisu, D. and Matsuda, Y. (2023+) Local polynomial trend regression for spatial data on $\mathbb{R}^d$. arXiv:2211.13467. Bernoulli in press.

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