Pythonで地理データを階層ベイズモデルで推定・可視化してみた

e-Stat APIをPythonから使って地理データを可視化してみた
e-Stat : https://www.e-stat.go.jp/SG1/estat/eStatTopPortal.doe-Stat API : http://www.e-stat.go.jp/api/e-Statは日本政府が調査した統計データを閲覧・ダウンロードできるよう管理されたポータルサイトです。このサイトからよく人口データなどをダウンロードし...

続き。

まぁ、続きと言いますか、少し前に書籍の『岩波データサイエンスVol.4 -地理空間情報処理-』を読みました。

そういえば、今は『Vol.5 -スパースモデリングと多変量データ解析-』も出ていますね。

で、Vol4の中で、e-Statから取得した地理データから階層ベイズモデルを用いて地域特徴を推定する事例を紹介されており、面白そうなのでやってみました。

ただ、上記書籍では、e-Statから普通にダウンロードしたデータを、Stan/BUGSなどのベイズ統計言語でデータを推定し、そこからまた別の地理データ可視化ソフトウェアを使って、推定結果を可視化していましたので、今回はそれらを全てPythonで一括してやってみようと思います。

問題と統計モデル

地域別の自殺リスクを推定する

上記の書籍では、地図上に表現された地理的な隣接情報を利用して、空間的な相関を考慮した階層ベイズモデルを用いた「地域別の自殺リスク」の推定を行って、その結果をコロプレスマップとして地図上に可視化、ということをやっています。

e-Statでダウンロードした地域別の自殺数を用いて、普通に自殺リスクを算出しようとすると、地域によっては値が小さすぎて、適切な自殺リスクが得られないため、ベイズモデルで空間的な相関を導入することで、これを解決するという方法を提案されています。

標準化死亡比(SMR)

「自殺リスク」とは何かという話ですが、これについては疫学の分野において、ある基準集団と比べて相対的にリスクが高いか低いかを表す指標で、標準化死亡比(SMR)というものがあります。

まず、ある地区が、基準集団と同じリスクを持つ場合の死亡数を、期待死亡数と呼び、下記で計算できます。

    \[y = \Sigma_i(\frac{Z_i}{N_i}{\times}{n_i})\]

Z_i : 基準集団での年齢iの死亡数
N_i : 基準集団の年齢iの人口
n_i : その地域の年齢iの人口

観測された死亡数が、この期待死亡数に比べて多い場合、死亡のリスクが大きいと考えられ、これを標準化死亡比(SMR)と言います。

    \[SMR=\frac{z}{y}\]

z : 観測死亡数

データからこの比率を単純に計算することもできますが、今回はこれを推定する方向で実装してみます。

統計モデル

上記を踏まえて、SMRを地域の相関が考慮されるように推定する階層ベイズモデルを次のように設計します。

    \[z_i{\sim}Poisson(\lambda_i)\]

    \[\lambda_i=y_i\exp(\phi_i+\psi_i)\]

    \[\phi_i{\sim}N(0, \sigma^2)\]

    \[\psi_i|\psi_{j{\neq}i}{\sim}N(\frac{1}{m_i}\Sigma_{j{\in}n(i)}\psi_i, \frac{1}{m_i}\sigma^2_{\psi})\]

z_i:地区iの死亡数
y_i:地区iの期待死亡数
\phi_i:地区i固有の効果を表すパラメータ
\psi_i:地区iは隣接地区と似通った傾向を持つことを表すパラメータ

データは死亡数なのでポアソン-対数正規モデルとし、ポアソン分布のパラメータに地域相関のあるSMRがかかる形にしています。

ここで、\exp(\phi_i+\psi_i)が地区iのSMRの推定値に相当します。

PyMC

PyMCは、ベイズ統計モデルの設定とMCMCサンプリングをPythonで実行できるライブラリです。

PyMC3 : https://github.com/pymc-devs/pymc3

PyMCにはPyMC2系列とPyMC3系列があるようです。

どっちがどうとかはよく分かりませんが、PyMC3の方が早いとかどこかで見ましたので、こちらを使うことにしました。

ちなみにPythonでのベイズ統計ライブラリといえば、Stanを実行できるPyStanもあります。

こちらでもよかったのですが、モデルの記述がまんまStanになりますので、ちょっと新鮮味を求めて今回はPyMCを使ってみました笑

比較してみると結構記述はPyMCの方がスマートで、直感的にモデルの構造が分かりやすい記述なのかなと思います。

ただ、インターネット上に情報がなく、PyStanは結局Stanですので、StanやRStanで書かれた情報とかも参考になる分、こちらは実装に苦労しました…。

PyMCの実装例

PyStanの実装例

地理データを階層ベイズモデルで推定・可視化

それでは、実装になります。

前回投稿同様に、データはe-Stat APIから直接取得します。

e-Stat API : http://www.e-stat.go.jp/api/

推定はPyMC3で行い、可視化は前回同様、Foliumを使って、そのまま分析環境で可視化してみます。

Folium : https://github.com/python-visualization/folium

コードが下記になります。

都道府県別でリスク推定をしてみた結果になります。

思ったよりも、地方の方がリスクが高いように推定されているみたいで、期待死亡数の方も、全体的に地方の方が高い傾向にありました。

ちなみに、書籍では市区町村別のデータにしていたため、おそらく欠損とかもある中、うまく推定をしているのでしょうが、そのデータがなぜか見つかりませんでした…。

本当は同じデータを使って、答え合わせしたかったのですが…。

それでも、実装の勉強になりました。

個人的には、地域相関を入れるCAR事前モデルの部分が、Stanですとパッケージで自動的にやってくれるみたいですが、PyMCはそれがないので自前で作らなければならず、そこが正しく出来ているのか不安です。

特に、多変量ガウス分布の分散共分散行列で地域相関を表現させる際の重みの調整はどうすれば良いか分からず、結局そこにも事前分布を入れて、回してみるようにしてみた、といったところです。

 

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です