はじめに
この記事はNTTテクノクロス Advent Calendar 2018の21日目の記事です．
こんにちは。NTTテクノクロスの奥野です。普段はhitoe®やAIを活用してバックオフィス業務を改革，レコメンドサービスShowBeeなどを中心にデータ分析業務に携わっています．
タイトルがかなり大げさですが，あまりクリギングについて書かれたブログを見かけないので宣伝がてら紹介してみます．
元ネタはOR学会機関紙（Vol.65, 2019年）に掲載（予定）されている新規店舗の売上予測の援助をGDPデータに適用してみたら予測精度がよくてクリギングって有用ですよ！って内容です．
クリギング
クリギングは地球統計科学でよく利用される手法で流行りの機械学習分野ではガウシアンプロセスと呼ばれています．もともとはD.G. Krigeによって提案された空間内挿手法の一つで距離の近い観測点のデータは大きな類似性を持つという空間相関構造を利用して，空間内挿を行う手法です．
クリギングにはデータの空間自己相関のみを考慮する通常クリギング (ordinary kriging) や線形回帰モデルの誤差の共分散を共分散関数でモデル化する普遍型クリギング（universal kriging）等がありますが，今回は普遍型クリギングの応用例について紹介します．
今回利用するデータについて
独立行政法人統計センターにて提供された教育用標準データセット（Standardized Statistical Data Set for Education: SSDSE）と内閣府が公開している県民経済計算（平成13年度 - 平成26年度）（1993SNA、平成17年基準計数）を利用させていただきました．
SSDSEはe-Statに収録されている「都道府県・市区町村のすがた（社会・人口統計体系）」の市区町村データから約100項目余を抜き出し、縦に市区町村、横にデータ項目が並ぶ表形式のデータに整備されたものです．
モデル化
今回は県庁所在地とそのGDPデータをクリギングで予測します．
クリギングは知ってるって方は「モデルの性能検証」までスキップしてください．
定式化1
地点 $s_{j} \in {j = 1, 2, \dots, n }$ における都道府県のGDPの対数を$y(s_{j})$とします．
このとき地点$s_j$は各都道府県の県庁所在地の位置を表します．
$y(s_{j})$が空間過程モデルで定義できると仮定し，$\boldsymbol{y}(\boldsymbol{s}) = (y(s_{1}), \ldots, y(s_{n}))^{\mathrm{T}}$を
\boldsymbol{y}(\boldsymbol{s}) \sim N(\boldsymbol{x} \boldsymbol{\beta}, \boldsymbol{\Sigma})
と定式化できます．ここで$\boldsymbol{x}$は$n \times k$の説明変数行列であり， $\boldsymbol{\beta}$ は$k \times 1$のパラメータベクトルを表します．
分散共分散行列は$\Sigma = \tau^{2}\boldsymbol{I} + \sigma^{2}\boldsymbol{H}(\phi)$と定式化します．
ここで$I$は $n \times n$ の単位行列，$\boldsymbol{H}(\phi)$はその$i, j$要素を相関関数で与える$n \times n$の相関行， $\tau^{2}, \sigma^{2}, \phi$ はそれぞれナゲット，
パーチャル・シル，レンジと呼ばれるパラメータです（詳細は参考文献をご参照してください）．
県庁所在地の距離が近いほど影響が高く，遠いほど影響が小さくなることを表すため$\Sigma$の${i, j}$要素である
$\Sigma_{i, j}$を指数型の共分散関数$C(d_{i, j})$を用いた場合，
C(d_{i, j}) = \exp(- \dfrac{d_{i, j}^{2}}{\phi^{2}})
で与えられます．ここで$d_{i, j}$は地点$s_{i}, s_{j}$のユーグリッド距離を表しており，等方性を仮定しています（異方的な場合は参考分を参照してください）．
そのため分散共分散行列は$C(d_{i, j}) = C(d_{j, i}) ^{\mathrm{T}}$と対称行列となります．
分散共分散のパラメータを推定するにあたり，パラメータ集合$\theta = (\beta^{\mathrm{T}}, \sigma^2, \tau^2, \phi)^{\mathrm{T}}$に関する事前分布は独立とします．ここで計算の都合上，式 $\boldsymbol{y}(\boldsymbol{s}) \sim N(\boldsymbol{x} \boldsymbol{\beta}, \boldsymbol{\Sigma})$ を次のように条件付きモデルに書き換えます．
\boldsymbol{y} \mid \boldsymbol{\beta},\boldsymbol{ \eta}, \tau^2 \sim N(\boldsymbol{x\beta} + \boldsymbol{\eta}, \tau^2 \boldsymbol{I}) \\
\boldsymbol{\eta} \mid \sigma^2, \phi \sim n(\boldsymbol{0}, \sigma^2 \boldsymbol{H}(\phi))
ここで，$\boldsymbol{\eta}$は空間ランダム効果$\eta(s_i)$からなる$n \times 1$のベクトルである．このように階層ベイズモデルとして定式化することで，
分散パラメータに関する条件付き時ご分布がMCMCによって計算することが可能となります．
事後確率密度関数は尤度関数と事前分布の積に比例するため
f(\boldsymbol{y} \mid \boldsymbol{\beta}, \boldsymbol{\eta}, \tau^2) f(\boldsymbol{\eta} \mid \sigma^2,\phi)p(\boldsymbol{\beta})p(\tau^2)p(\sigma^2)p(\phi)
と定式化できます．
目的変数と説明変数について
地点$s_i$を県庁所在地とし，目的変数を$y_i$を2014年度の県内総生産とします（2014年を選択した理由については後述）．
$\boldsymbol{x}$である説明変数にはSSDSEのデータを利用しました．説明変数として利用する具体的な項目は表に示すとおりです．
データは項目によって毎年・隔年・5年周期取得されているためすべての項目が同年度で取得されているわけではないため，
データ項目数が一番多い2014年度のデータを利用しました．
※ モデル選択として複数の組み合わせからAICが最小となるものを選んでいます．
| 記号 | 説明変数 | 取得方法 |
|---|---|---|
| $y_{s_{j}}$ | 県内装生産 | 県民経済計算（平成13年度 - 平成26年度）（1993SNA、平成17年基準計数） |
| $x_{s_{j}, 1}$ | 事業所数[10万所] | SSDSE |
| $x_{s_{j}, 2}$ | 従業員数[10万人] | SSDSE |
| $x_{s_{j}, 3}$ | 小売店舗数[万店舗] | SSDSE |
| $x_{s_{j}, 4}$ | 飲食店舗数[万店舗] | SSDSE |
| $x_{s_{j}, 5}$ | 大型小売店数[万店舗] | SSDSE |
| $x_{s_{j}, 6}$ | 医師数[万人] | SSDSE |
| $x_{s_{j}, 7}$ | 歯科医師数[万人] | SSDSE |
| $x_{s_{j}, 8}$ | 薬剤師数[万人] | SSDSE |
推定結果とその活用方法2
パラメータの推定結果
以上のデータから推定したパラメータを表に示します．
表中には事後平均と標準偏差の目安として2.5%と97.5%の領域を示しておきました．
表を見てみると「小売店舗数」のみが95\%信頼区間の意味で有意であることがわかります．
ただし共分散関数のパラメータの推定値から誤差の17.6\%($=\tau^{2} / (\sigma^{2} + \tau^{2})$)は空間相関成分であることから推定に利用した説明変数では説明できていない要因（例えば輸入量，輸出価格，為替レート等）が多く存在することを示唆しています．
GDPの変動要因として民需・公需・外需があげられるが，本モデルでは民需が取り込まれており，その要因である小売店舗数が有意な結果となったと考えられます．
| 説明変数 | 2.5% | 50.0% | 97.5% |
|---|---|---|---|
| 切片 | 14.06 | 14.41 | 14.77 |
| 事業所数 | -2.86 | 0.02 | 3.02 |
| 従業員数 | -0.07 | 0.05 | 0.16 |
| 小売店舗数 | 0.00 | 1.13 | 2.25 |
| 飲食店舗数 | -2.63 | -1.24 | 0.12 |
| 大型小売店数 | -17.33 | 9.71 | 37.84 |
| 医師数 | -1.64 | 0.35 | 2.39 |
| 歯科医師数 | -7.65 | -2.67 | 2.42 |
| 薬剤師数 | -1.83 | -0.24 | 1.36 |
| $\sigma^{2}$ | 0.09 | 0.14 | 0.24 |
| $\tau^{2}$ | 0.01 | 0.03 | 0.07 |
| $\phi$ | 3.56 | 17.78 | 28.67 |
モデルの性能検証
推定結果から小売店舗数が有意な結果となったと示しましたが，モデルが妥当であるかはわからない状態です．
そこでモデルの妥当性として今回は予測性能の良さで検証したいと思います．
性能の検証方法は県月GDPを47分割交差検定(47-fold cross-validation)で検証します．
通常の47分割交差検定法は，データ全体を無造作に47組に分割し，『46組のデータを用いてパラメータ推定を行い，残り1組のデータを予測して精度検証を行う』という過程を47回繰り返し予測精度を評価します．
また精度検証にあたり比較対象手法として通常クリギング
\boldsymbol{y}(\boldsymbol{s}) \sim N(\boldsymbol{0}, \boldsymbol{\Sigma})
を利用しました．通常クリギングの推定も同様にMCMCにて 15,000 回の サンプリングを行い，
はじめの 10,000 回をバーンイン期間として破棄し，残りの 5,000 個の MCMC 標本に基づいて得た事後統計量の平均値を推定値としています.
47 分割交差検定によるGDPの予測精度を表に示す．表中に示す値は$\ln$(予測したGDP)と
$\ln$(GDP) の平均二条平方根誤差（RMSE）を計算し，$\exp$(RMSE)と変換した値を示します．
$\exp(RMSE) = 1.0$ なら予測したGDPが完全に一致することを表し，
$\exp(RMSE) = 1.1$ なら予測したGDPが平均的に10%異なる状態を示します．
| 手法 | RMSE | exp(RMSE) |
|---|---|---|
| 通常クリギング | 0.72 | 2.07 |
| 普遍型クリギング | 0.069 | 1.07 |
表を見ていただければ結果は歴然ですね．通常クリギングと比べて高い精度で予測することが可能となっています．
オーバーフィットの問題はおいておいて推定誤差も10%以下であることからモデルの予測性能が高く実務で適用できる可能性が期待されます．
おわりに
都道府県別GDPを利用して情報を抽出する方法論として空間相関を構造化して
内挿を実行できる普遍型クリギングに着目し，SSDSEに対して実証実験を行い，
普遍クリギングのGDP予測（内挿）の適用可能性を検証しました．
その結果，県庁所在地において10%以下の高精度な予測が可能であることが明らかになり，
空間的な予測が可能であることを確認した．またGDPに寄与する要因の影響度をモデルの大域変動項である
説明変数の回帰係数を推定することで明らかにしたということでタイトルの「予測と情報抽出」を少しは達成できたかな？
この見づらい数式のスクロールバーを消す方法がなにかあるのだろうか．．．