16. OpenBUGSで欠測値補完

藤宮

『はじめに』

今回は、OpenBUGS ( http://www.openbugs.net/ ) を用いた欠測値補完の例です。OpenBUGSは、現在でもバージョンアップが進められています。また、RからBRugsパッケージを利用することで、完全にバックグラウンドで処理させることができます。WinBUGSだとRから起動されたときにいったんマウスなどの制御がWinBUGSのウインドウに取られてしまう問題があったりしますが、OpenBUGSであれば、MCMCの処理中にフォアグランドでメールを見たり作成したりしていても問題なく気にせず使えるので便利です。(PCによっては遅くなりますが)。データとしては、どなたでもダウンロードして利用できる人口推計データです(統計局: http://www.stat.go.jp/ )。ここでは、人為的に欠測値を作り、推定結果と比較してみます。年度別、年齢別の2次元データです。このようなデータで、任意の場所の欠測値をしっかり補完するためには、条件付きの循環参照を実現しなければなりません。したがって残念ながら現在のJAGSは使えません(ver. 3.4.0時点において)。

 

『シミュレーション用人口推計データ』

図1に統計局から取得した人口推計データを示します。0まで切れ込んでいる部分は意図的に欠測値とした部分です。26歳~30歳までの5年間と、60歳以上ぐらいに5箇所程度欠測値にしてみました。

tech016_01

図1 欠測値を持つ人口データ(0になる深い切れ込みが欠測値)

 

『MCMCのモデル』

人口推計データは、1年毎に全員が1歳年を取りますので、2005年から2006年に移るときに生存率分だけの人数を斜め方向にシフトさせた関係になります。逆に見れば、生存率で割った数だけ前の年の1歳若い側の人口に戻せるという関係です。どこにでも欠測値を許すため、両側からその関係を記述します。さらに生存率関数の事前分布として、下記のcloglog関数を使用しました。生存率としているため、多少変形しています。この生存率のハイパーパラメータをサンプリングしながら、読み込んだデータの欠測値も更新しつつ最大尤度になるようにMCMCが進みます。事前分布である生存率のハイパーパラメータはロケーション、スケールの2個があります。

tech016_02

図2 生存率の事前分布関数の例(1-icloglog関数)
(徐々に生存率が落ちながら右下では急激に0になる)

 

『結果』

MCMCを実行した結果を図3以降に示します。図4はサンプリング系列です。推定している欠測値や生存率の事前分布関数に使用したハイパーパラメータ2個(beta0, beta1)も非常に安定したサンプリング結果が得られました。欠測値を推定した結果の相対誤差を図5に示します。25~30歳においてまとめて欠落させた部分の相対誤差が2~4%です。その他に100歳付近で人数が少なくなるため、相対誤差が大きくなってしまいました。

tech016_03

図3 欠測値補完の結果

 

tech016_04

tech016_05

図4 生存率曲線のハイパーパラメータのサンプリング状況(3回分を重ねてプロット)

 

tech016_06

図5 オリジナルのデータに対する推定結果の相対誤差