空間データのためのブートストラップ法

東京工業大学/株式会社Nospare の栗栖です．
前の記事 (「空間データのためのブートストラップ法」(解説編)) では空間データ・時空間データ分析に関する最近の研究成果 (Kurisu, Kato and Shao(2021)) で提案した空間ワイルドブートストラップ(spatially dependent wild bootstrap, SDWB)について解説しました．この記事では (解説編) の続きとして数値実験と降水量データの分析例について紹介します．紹介する分析は全て R を利用して行っています．

数値実験

データ生成過程

以下のモデルから数値実験のデータを生成します：

\boldsymbol{Y}(\boldsymbol{s}_{i}) = A\boldsymbol{F}(\boldsymbol{s}_{i}) + \boldsymbol{R}(\boldsymbol{s}_{i}),\ i=1,\cdots,n=100,

$A = (A_{j,k})$：$p \times 5$ 行列，$A_{j,k}$ は独立に $[-1,1]$ 上の一様分布から生成，
$\boldsymbol{F}$：各成分が独立な平均 $\boldsymbol{0}$ の $5$-変量ガウス過程，$\boldsymbol{F}(\boldsymbol{s}) = (F_{1}(\boldsymbol{s}),\cdots,F_{5}(\boldsymbol{s}))'$，$\boldsymbol{s} \in \mathbb{R}^{2}$．
特に $F_{j}$ は全て以下の共分散関数をもつとします¹：

\text{Cov}(F_{j}(\boldsymbol{s}), F_{j}(\boldsymbol{0})) = (3\|\boldsymbol{s}\|+1)e^{-3\|\boldsymbol{s}\|}, \|\boldsymbol{s}\|=\sqrt{s_{1}^{2} + s_{2}^{2}}

$\boldsymbol{R}(\boldsymbol{s}_{i})$：$p$-変量の標準正規分布に従う独立同分布な観測誤差．
$\boldsymbol{s}_{i} \in R_{n} = [-\lambda_{n}/2, \lambda_{n}/2]^{2}$，$\boldsymbol{s}_{i}$ は独立に $R_{n}$ 上の一様分布から生成．
$p = 10,100,400$，$\lambda_{n} = 15,25$．

$n=100$ なので $p=100, 400$ が高次元の場合に対応．
空間ワイルドブートストラップの実装に必要になるガウス過程のカーネル関数としては

a(x) = 
\begin{cases}
1-|x| & \text{$|x|\leq 1$ のとき} \\
0 & \text{それ以外}
\end{cases}

さらに $b = 1,2,\cdots ,10$ とします．

空間ワイルドブートストラップの計算

解説編でも書きましたが，空間ワイルドブートストラップの計算方法についてもう一度確認しておきます．

(Step1) 標本平均 $\overline{\boldsymbol{Y}}_{n} = {1 \over n}\sum_{i=1}^{n}\boldsymbol{Y}(\boldsymbol{s}_{i})$ と $\hat{\Sigma}^{\boldsymbol{V}_{n}}$ を計算：

\hat{\Sigma}^{\boldsymbol{V}_{n}} = {1 \over n^{2}\lambda_{n}^{-d}}\sum_{i_{1}, i_{2} = 1}^{n}(\boldsymbol{Y}(\boldsymbol{s}_{i_{1}}) - \overline{\boldsymbol{Y}}_{n})(\boldsymbol{Y}(\boldsymbol{s}_{i_{2}}) - \overline{\boldsymbol{Y}}_{n})'a\left({\boldsymbol{s}_{1} - \boldsymbol{s}_{2} \over b_{n}}\right) 

(Step2) ブートストラップ標本の数 $B$ を設定し，共分散行列 $A = (a(\boldsymbol{s}_{i}-\boldsymbol{s}_{k})/b_{n}))_{1 \leq i,k \leq n}$ をもつ平均 $\boldsymbol{0}$ の $n$ 次元正規分布に従う変数 $\boldsymbol{W}^{(1)}$$,\dots,$ $\boldsymbol{W}^{(B)}$ を生成：

\boldsymbol{W}^{(k)} = (W_{1}^{(k)},\dots,W_{n}^{(k)})',\ k=1,\dots, B.

(Step3) ブートストラップ標本 $\mathcal{Y}^{(1)}_{n}, \dots, \mathcal{Y}^{(B)}_{n}$ を生成：

\begin{align}
\mathcal{Y}^{(k)}_{n} &= \{\boldsymbol{Y}^{(k,\ast)}(\boldsymbol{s}_{1}),\dots, \boldsymbol{Y}^{(k,\ast)}(\boldsymbol{s}_{n})\},\\
\boldsymbol{Y}^{(k,\ast)}(\boldsymbol{s}_{i}) &= \overline{\boldsymbol{Y}}_{n} + \left(\boldsymbol{Y}(\boldsymbol{s}_{i} ) - \overline{\boldsymbol{Y}}_{n}\right)W^{(k)}_{i},\ k=1,\dots,B.
\end{align}

(Step4) ブートストラップ標本 $\mathcal{Y}_{n}^{(k)}$ ごとに

T_{n}^{(k)} = \sqrt{\lambda_{n}^{d}}\max_{1 \leq j \leq p}\left|{\overline{Y}_{j,n}^{(k,\ast)} - \overline{Y}_{j,n} \over \sqrt{\hat{\Sigma}_{j,j}^{\boldsymbol{V}_{n}}}}\right|,\ k=1,\dots, B

を計算して $T_{n}^{(1)},\dots, T_{n}^{(B)}$ の経験分布から分位点 $q_{n}(1-\tau)$ の推定値 $\hat{q}_{n}(1-\tau)$ を計算．

数値実験の結果

以下の図は上記のデータ生成過程から得られたデータをもとに空間ワイルドブートストラップを利用して 95% 信頼区間を1000回繰り返し計算し，

{\text{1000回中 $p$ 個の信頼区間がすべて $0$ を含んでいた回数} \over 1000}

(経験被覆確率，empirical coverage probability) を計算したものです．この例ではブートストラップ標本の数 $B = 1500$ としました．今回は95%信頼区間を計算してるので，図の中の折れ線が黒の横線 (coverage probability = 0.95) に近いほどパフォーマンスが良いことを示しています．

バンド幅 $b$ の選択について：
図を見ると最適な $b$ が1個あるというよりは $b$ を適当な範囲 (上の例だと$3 \leq b \leq 7$ くらい) の中で選んでおけばそれなりに満足のいくパフォーマンスが保証される (折れ線の値が 0.95 に近い) ことがわかります．このことから，実際のデータ分析では
(1) 観測地点間の距離を計算するために $\lambda_{n}$ をデータの分析者が好きなように決めて，
(2) 複数の $b$ に対して分析を行い，それぞれの $b$ ごとの結果を比較し，
(3) 結果があまり変化せず安定している $b$ の範囲から $b$ を1つ選べばよい
ということを示唆しています．
また数値実験の結果から，$\lambda_{n}=15$ のときより $\lambda_{n}=25$ の時の方がパフォーマンスが良くなっていることもわかります．一方（解説編）で述べたように，標本平均 $\overline{\boldsymbol{Y}}_{n}$ が $\boldsymbol{\mu}=E[\boldsymbol{Y}(\boldsymbol{s})]$ を近似する精度が $\lambda_{n}^{-d/2}(\log p)^{1/2}$，(即ち，$\overline{\boldsymbol{Y}}_{n}-\boldsymbol{\mu}=O_{p}(\lambda_{n}^{-d/2}(\log p)^{1/2})$ であることがわかります．これは $\lambda_{n}$ が大きい時の方が理論的にもパフォーマンスが良くなることを意味しますが，数値実験でもこのことが確認できたことになります．

データ分析

以下では空間ワイルドブートストラップをアメリカの降水量のデータの分析に適用した結果について紹介します．

##アメリカ降水量データ

下の図はアメリカの中西部の州(ノースダコタ，サウスダコタ，ネブラスカ，ミネソタ，ウィスコンシン，アイオワ) において，1919年-1994年の間の各年の年間降水量の観測値が利用可能な気象台の配置図です (この例では観測地点の数 $n=101$)．ここでは各地点 $\boldsymbol{s}$ の年間降水量 $M_{1919}(\boldsymbol{s}),\cdots,M_{1994}(\boldsymbol{s})$ がどのように変動しているのか分析を行ってみたいと思います²．具体的には各地点 $\boldsymbol{s}$ で観測された 76年間 (1919年-1994年) の年間降水量を $Y(t,\boldsymbol{s})=\log(1+M_{t}(\boldsymbol{s}))$, $t=1919,\dots,1994$ と変換した値の差 $Y_{1}(\boldsymbol{s}) = Y(1920,\boldsymbol{s}) - Y(1919,\boldsymbol{s})$，$\cdots$，$Y_{75}(\boldsymbol{s}) = Y(1994,\boldsymbol{s}) - Y(1993,\boldsymbol{s})$ を $75$ 次元の空間データ $\boldsymbol{Y}(\boldsymbol{s}) = (Y_{1}(\boldsymbol{s}),\cdots,Y_{75}(\boldsymbol{s}))'$ に変換し，$75$ 個の仮説検定

\mathbb{H}_{0}: \mu_{j} = E[Y_{j}(\boldsymbol{s})] = 0\quad \text{v.s.}\quad \mathbb{H}_{1}: \mu_{j} = E[Y_{j}(\boldsymbol{s})] \neq 0

を実行してみます．これは各年で前年と比べ降水量が変化したかどうかを同時に検定していることに対応しています．

年間降水量の分析

今回の例では $\lambda_{n}=15$ とし，空間ワイルドブートストラップの実装に必要になるガウス過程のカーネル関数として数値実験と同じ

a(x) = 
\begin{cases}
1-|x| & \text{$|x|\leq 1$ のとき} \\
0 & \text{それ以外}
\end{cases}

(Bartlett カーネル) を使います．さらに $b = 1,2,\cdots ,10$ として分析してみます．下の図はそれぞれ上から順に $b=1,4,7,10$ の場合における各年における同時95% (濃い灰色)・99% (薄い灰色)信頼区間を描いたものです．この例でもブートストラップ標本の数 $B = 1500$ としました．
全ての図において黒の実線が各年における $\mu_{j}$ の推定値で，仮説を棄却するかどうかの基準となる線を青色で描いています．この青い線から信頼区間が上または下に外れれば帰無仮説は棄却されます．図を見ると $b$ の値によって信頼区間の幅が変化していることがわかります．実際 $b=1 \to 4 \to 7 \to 10$ の順に棄却された帰無仮説の数は $39 \to 24 \to 21 \to 24$ と変化します.

さらに下の図は $b=1,2,\cdots,10$ の全てにおいて有意水準99%で棄却された帰無仮説を集め，対応する年をプロットしたものです．各年で99%信頼区間がプラスに外れた年 (信頼区間の下限がプラスの年) は "up"，マイナスに外れた年 (信頼区間の上限がマイナスの年) は "down" としてプロットしています．"up" または "down" となった年は75年という長期間でみても年間降水量の変動が大きかった年であることを意味します．以上の結果により，今回利用したデータが得られた地域では年間降水量は時間方向に非定常性をもつことが示唆されます．

論文の中では他の数値実験の結果についても報告していますので興味のある方はご参照ください．

まとめ

この記事では Kurisu, Kato and Shao(2021) で提案した空間ワイルドブートストラップを利用した時空間データ分析の例について紹介しました．今後は他の時空間統計の話題についても記事を書いていこうと思います．
株式会社Nospareには空間統計・時空間統計に限らず，統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては 株式会社Nospare までお問い合わせください．

参考文献
[1] Kurisu, D., Kato, K. and Shao, X. (2021) Gaussian approximation and spatially dependent wild bootstrap for high-dimensional spatial data. R&R at Journal of the American Statistical Association. arXiv:2103.10720.