はじめに
東京大学/株式会社Nospare リサーチャーの栗栖です．
この記事は最近の研究成果「空間データに対する局所多項式トレンド回帰」(Kurisu and Matsuda(2024))についての解説記事の続編です．前回の記事では本研究の目的と局所多項式回帰推定量 (local polynomial (LP) estimator) の漸近分布について紹介しましたが，この記事ではLP推定量の信頼区間の構成方法について紹介します．
モデル
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(x): 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$ をLP推定量で推定し，信頼区間を構成する問題を考えます．推定にあたり，上記のモデルに対して以下を仮定します．
$m(x)$ の LP推定量は以下で与えられることが分かります．
\begin{align*}
\hat{\beta}(x) &=\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 \\
&= (\hat{m}(x),\partial_1 \hat{m}(x),\partial_2 \hat{m}(x),\dots)'
\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}
です．詳細については前回の記事を参照してください．
LP推定量の漸近正規性
次にLP推定量の漸近的性質について紹介しておきます．適当な条件の下で，$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*}
ここで，$\sigma_e(x)=E[e(0)e(x)]$ は空間相関をもつ誤差項の相関関数，$f(x)$ は $x$ における $\boldsymbol{X}_i$ の密度関数の値であり， $S$，$V$ はカーネル関数 $K$ に依存する行列です．また
\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*}
です．これらの詳細な定義については前回の記事を参照してください．
LP推定量の一様収束レート
LP推定量については適当な集合上での一様な収束を示すことができます．特に $m(x)$ に対して，適当な条件のもとで以下が成り立ちます：
\begin{align*}
\sup_{x \in [-1+h/2, 1+h/2]^2}|\hat{m}(x) - m(x)| &= O_p\left(h^{p+1} + \sqrt{\log n \over Ah^2}\right).
\end{align*}
上記の結果以外にも，$m(x)$ の偏微分の値についても一様収束レートを導出することができます．
LP推定量の漸近分散の推定
上記の結果において，推定量の漸近分布として現れる正規分布の分散を推定する問題を考えましょう．漸近分散の推定量として，以下の推定量を考えます：
\begin{align*}
\hat{W}_{n}(x) &= \left(\int K^2(u)du\right)^{-1}{\hat{W}_{n,1}(x) \over \hat{g}^2(x)}.
\end{align*}
ここで，
\begin{align*}
\hat{W}_{n,1}(x) &= {A \over n^2 h^2}\sum_{i,j=1}^{n}K\left({\boldsymbol{X}_i/A - x\over h}\right)K\left({\boldsymbol{X}_j/A - x \over h}\right)\\
&\quad \times \left(1-{\|\boldsymbol{X}_i - \boldsymbol{X}_j\| \over b}\right)1_{[0,1]}\left({\|\boldsymbol{X}_i - \boldsymbol{X}_j\| \over b}\right)\\
&\quad \times (Y_i - \hat{m}(\boldsymbol{X}_i/A))(Y_j - \hat{m}(\boldsymbol{X}_j/A)),\\
\hat{g}(x) &= {1 \over nh^2}\sum_{i=1}^{n}K\left({\boldsymbol{X}_i/A - x \over h}\right)
\end{align*}
であり，$b$は $n \to \infty$ として $b \to \infty$ となる定数列とします (ただし $b < Ah$)．
上記の推定量 $\hat{W}_n(x)$ は時系列データ分析において，HAC推定量 (heteroskedasiticity and autocorrelation consistent estimator) (データの自己相関を考慮した統計量の漸近分散の推定量) として知られている推定量の構成のアイデアを用いています．
このとき，適当な条件のもとで $n \to \infty$ として以下が成り立ちます：
\begin{align*}
\hat{W}_n(x) &\stackrel{p}{\to} {2\kappa \over f(x)} + \int_{\mathbb{R}^2}\sigma_e(x)dx.
\end{align*}
この結果を利用すれば，例えば$m(x)$ の $95$%信頼区間を以下のように構成することができます：
\begin{align*}
C_{0.95}(x) =\left[\hat{m}(x) - 1.96\sqrt{\hat{W}_n(x)[S^{-1} V S^{-1}]_{1,1} \over Ah^2}, \hat{m}(x) + 1.96\sqrt{\hat{W}_n(x)[S^{-1} V S^{-1}]_{1,1} \over Ah^2}\right]
\end{align*}
ここで，$[S^{-1} V S^{-1}]_{1,1}$ は $S^{-1} V S^{-1}$ の $(1,1)$成分です．実際以下を示すことができます：
\begin{align*}
\lim_{n \to \infty}P\left(m(x) \in C_{0.95}(x)\right) = 0.95.
\end{align*}
回帰関数の信頼曲面の構成
ここまでは回帰関数の各点での信頼区間を構成する方法について説明しました．一方で，回帰関数は2次元空間内の集合$R=[-A,A]^2$上で定義された関数なので，$R$上で回帰関数がどのような形状をしているのかにも興味があります．$R$上で回帰関数の形状とその推定精度を評価する方法としては以下の方法が考えられます．
[1] $R$内の有限個の$L$個の代表点$\{x_\ell\}_{\ell=1}^L$をとる．
[2] 各点$x_\ell$で信頼区間を構成
\begin{align*}
C_{0.95}(x_\ell) =\left[\hat{m}(x_\ell) - \xi_{0.95}\sqrt{\hat{W}_n(x_\ell)[S^{-1} V S^{-1}]_{1,1} \over Ah^2}, \hat{m}(x_\ell) + \xi_{0.95}\sqrt{\hat{W}_n(x_\ell)[S^{-1} V S^{-1}]_{1,1} \over Ah^2}\right]
\end{align*}
ここで，$\xi_{0.95}$は$L$個の独立な標準正規分布$Z_1,\dots,Z_L$の絶対値の最大値の分布の0.95分位点で，特に$\xi_{0.95}$は以下を満たす値です．
P\left(\max_{1 \leq \ell \leq L}|Z_\ell|>\xi_{0.95}\right) = 0.05.
[3] $\{C_{0.95}(x_\ell)\}_{\ell=1}^L$を線形補完して信頼曲面を描く．
以上の手続きにより回帰関数(空間トレンド)$m$の95%信頼曲面を構成することができます．特に，適当な条件の下で以下の結果を示すことができます．
\lim_{n \to \infty}P\left(m(x_\ell) \in C_{0.95}(x_\ell), \ell=1,\dots,L\right) = 0.95.
下の図は信頼曲面のイメージ図です．
次回の記事では数値実験の結果を通して理論的な結果が実際に成立していることを確認します．
まとめ
この記事では空間データに対するノンパラメトリック回帰分析に関して，最近の研究成果である Kurisu and Matsuda (2024) を紹介しました．株式会社Nospareにはノンパラメトリック分析や空間統計に限らず，統計学の様々な分野を専門とする研究者が所属しています．統計アドバイザリーやビジネスデータ分析につきましては株式会社Nospareまでお問い合わせください．
参考文献
[1] Kurisu, D. and Matsuda, Y. (2024) Local polynomial trend regression for spatial data on $\mathbb{R}^d$. arXiv:2211.13467. Bernoulli in press.