新規作成:2015年12月13日
最終更新:2016年9月22日

ローカルレベルモデルくらいでしたらdlmパッケージでも簡単に推定できます。
せっかくStanを使うのですから、もっと複雑なモデルを作ってみましょう。
今回は回帰係数が時間によって変化するモデルを作成します。

なお、今回の記事はちょっと複雑なモデルを作成する都合上、実験的な要素もあります。
そのため、ほかの記事に比べると、記載内容の精度が落ちております。この記事を参考にされる場合は、この点、注意してください。

すべてのコードをまとめたものはこちらにおいておきました。コピペする際はこちらを使用してください。

この記事はベイズ推定を応用して状態空間モデルを推定する一連の記事の一つです。
記事の一覧とそのリンクは以下の通りです。
ベイズ統計学基礎
ベイズと統計モデルの関係
ベイズとMCMCと統計モデルの関係
Stanによるベイズ推定の基礎
Stanで推定するローカルレベルモデル



スポンサードリンク

 

目次

1.今回扱うモデルの説明
2.シミュレーションデータの作成
3.Stanによる状態空間モデルのベイズ推定
4.結果の答え合わせ
5.推定結果の図示

 

1.今回扱うモデルの説明

Stanによる状態空間モデル推定の応用編として、今回は“慣れ”をモデリングします。

最初は効果があったのに、時間がたつにつれて効果が薄くなってくること、よくありますよね。
その「時間がたつにつれて効果が薄くなる」ことを統計モデルで表します。
具体的には、「回帰係数が時間によって変化する」統計モデルを作ります。

例えば、血圧を下げる薬があったとします。
血圧 ~ 切片 + 回帰係数×薬
上記の式は普通に回帰分析をすれば推定できるのですが、この回帰係数を時間によって変化するようにしてやりましょう。
そうすれば「徐々に回帰係数が0に近づく→薬に慣れて効果が薄くなる」ことを表すことができます。

ローカルレベルモデルだけですとなんだかベイズ推定のありがたみがありませんでした。
Stanは本来、今回のモデルのように複雑なモデルを推定するのにつかわれるべきでしょう。
Stanを使うことによって、現象を非常に柔軟にモデリングすることができます。
Stanによるベイズ推定の優秀性を確認いただければと思います。
そして、複雑なモデルであっても、基本はローカルレベルモデルと同じだということも感じ取ってもらえれば幸いです。

 

2.シミュレーションデータの作成

今回は血圧を下げる薬の効果が徐々になくなる(薬に慣れてしまう)ことをモデル化します。
しかし、血圧の実データは持っていないですので、シミュレーションして作ります。
結果の見やすさのため、少々現実離れした設定にしてあります。ご容赦ください。
また、ややこしくなるのを防ぐため、収縮期血圧だけを対象としました。

ちなみに、期間は110日間。
最初の10日は薬を入れていません。
後の100日に薬を投与したという前提で進めていきます。

まずは、薬を一切投与しなかった時の血圧の値をシミュレーションしました。
ここは単純にランダムウォークを仮定しています。
この値はdlmパッケージに合わせて「mu」という名前にしています。

血圧の「状態」の初期値を160として、そこから、平均0、標準偏差3の正規分布に従う正規乱数を足し合わせています。
ランダムウォークがわからない方はこちらの記事をご覧ください。

シミュレーションの結果はこちら。
データを表示させた後に描画しました。

時変係数モデル_1_血圧の状態
(画像をクリックすると大きくなります)

次はこの「ランダムウォークしている血圧 mu 」に薬の影響を入れていきます。
以下のようなモデルになります。
血圧 ~ 切片(ランダムウォークしている) + 傾き×薬の量
薬を使うと、血圧が下がることを想定しているので、傾きはマイナスになるはずです。

今回はさらに傾き(回帰係数)が時間によって変化するようにします。
”慣れ”をモデル化するという名目なので、薬の回帰係数が徐々に0に近づくようにしてやりましょう。
そこで、回帰係数にローカル線形トレンドモデルを用います。
ローカル線形トレンドモデルについてわからない方はこちらを参照してください。
こちらでは簡単な説明をするにとどめます。

ローカル線形トレンドモデルは、以下のように考えます。
1.薬の効果には、毎回「トレンド」の値が追加される
例:トレンドが常に「+1」のときは、薬を入れても毎日1ずつ「血圧が下がりにくくなる」ことになる
10日後には、「血圧が10下がりにくくなる」ことになる
初期に「血圧を―25する効果」があったとしたら、10日後には「血圧を―15する効果」でしかなくなっているというわけ。
これが”慣れ”ですね。
2.「薬の効果」はランダムウォークしている
3.「トレンド」の値もランダムウォークする
ランダムウォークするのは、「薬の効果」と、「薬の効果のトレンド」の2つであることに注意してください。

まずは、「薬の効果のトレンド」からシミュレーションします。

「薬の効果のトレンド」は初期値が0.005で、そこから平均0で標準偏差0.03の正規分布に従う乱数が足しあわされていったものだと仮定しています。
どんな結果になったのか確認してみます。

時変係数モデル_2_薬の効果のトレンド
(画像をクリックすると大きくなります)

トレンドの合計値は18です。
100日後には、「血圧が18だけ”下がりにくくなる”」ことがわかります。
本当はこんなあっけなく薬が効かなくなることはないと思います。
結果の見やすさのため、少々現実と違う前提にしてあることはご容赦ください。

次は、薬の効果のシミュレーションです。
先ほどシミュレーションされた「薬の効果のトレンド」が足しあわされていることに注意してください。

図示します。

時変係数モデル_3_薬の効果
(画像をクリックすると大きくなります)

最初は、薬を入れると、血圧が25も小さくなっていました。
しかし、徐々に慣れていくので、100日後には、血圧を7~8くらいしか下げる効果がないようです。

実際の薬の効果は、さらにノイズが加わると仮定しましょう。
最終的な薬の効果の時間的な変化は、以下のようになります。

図示します。

時変係数モデル_4_ノイズ付きの薬の効果
(画像をクリックすると大きくなります)

黒い線が、実際の薬の効果。赤い線がノイズをの取り除いた薬の効果を表しています。「トレンド」という言葉が少々紛らわしいのですが、ノイズや周期変動を取り除いた結果のこともトレンドといいます。
これでようやく「時間によって変化する薬の効果」をシミュレーションできました。
あとは、昔作った「ランダムウォークする血圧 mu」に薬の効果を足してやるだけです。

最初の10日は薬なしで、11~70日は薬あり。
効きが悪くなったので、71日後に薬の量を倍にし、101日後には初期の3倍の量にまで増やしたことにします。

どんな結果になったか、確認してみます。
薬の量を変えた日には縦線を入れておきました。

時変係数モデル_5_血圧プロット1
(画像をクリックすると大きくなります)

なんとなくそれっぽい結果になっていますね。
最初は薬を入れると、血圧がガクッと下がるんですが、徐々に血圧が上がっていきます。
薬を増やすと血圧は下がるんですが、徐々に効果がなくなっていく。
実際はここまで顕著ではないでしょうか、勉強のためのシミュレーションデータとしてはまずまずです。

せっかくなので、「状態」の値も重ね合わせてみましょう。

時変係数モデル_5_血圧プロット2
(画像をクリックすると大きくなります)

この緑線だの赤線だのの「状態」を今からStanを使って推定するわけです。
やってみましょう。

 

3.Stanによる状態空間モデルのベイズ推定

Stanでベイズ推定します。
やり方は以前と同じですので、ごく簡単に解説します。

ステップ1:Stanを使用する準備

まずは、Stanを使用する準備からです。
Stanをインストールした後に、以下を実行してください。
Stanのインストール方法などは以前の記事を参考にしてください。

ステップ2:統計モデルの記述

統計モデルに関しては、ローカルレベルモデルよりも大分と複雑になっています。
ですが、今回に限ってはそれほど難しくありません。
なぜならば、シミュレーションデータの作成とまったく同じ式を書くだけだからです。

モデルの記述は以下のようになります。

dataブロック

dataブロックでは、サンプルサイズや血圧の値、薬の量を入れています。

parametersブロック

parametersブロックはちょっと長いんですが、ローカルレベルモデルで推定されたものと解釈は同じです。
ローカル線形トレンドモデルや時変係数モデルは、ローカルレベルモデルを足し合わせたものなので。

ちなみに、データは自分で作りましたので、正しいパラメタはわかっています。
一覧を記していおきます。
muZero:160
sigmaV:10(標準偏差です。シミュレーションの最後の段階でノイズを加えるときにsd=10と出てきています)
sigmaW:3(標準偏差です)
coefMedicineTrendZero:0.005
coefMedicineZero:-25
sigmaWcoef:0.5(標準偏差です)
sigmaWcoefTrend:0.03(標準偏差です)
sigmaVcoef:2(標準偏差です)

modelブロック

大変長くて恐縮なのですが、シミュレーションデータ作成時とほぼほぼ同じfor分が出てきているので、対応はつかめるかと思います。
ちなみに、標準偏差の値が軒並み100で割られています。これは、標準偏差がとても小さな値になることがあるからです(例えばsigmaWcoefTrendは0.03です)。
あんまりにも絶対値が小さい値を推定するのはちょっと難儀なので、補正してみました。
ですので、標準偏差に関しては、推定された値を100で割ったものが正確な標準偏差になることには注意してください。

ステップ3:データの指定

Stanに渡すデータを作ります。
結果はこんな感じ。

ステップ4:ベイズ推定の実行

ベイズ推定します。
ちなみに、この計算は、管理人のPCで2.5時間ほどかかりました。
お出かけ前や、寝る前に計算されることをお勧めします。

コメントにも書いたのですが、警告メッセージが出てしまいました。
警告が出ないようにとcontrolを追加したのですが、ちょっとうまくいかなかったのでそのまま載せています。
本当はあまりよくないのですが、今回は答えがわかっているデータを相手にして、答え合わせができますので、そのままにしておいています。
今後、警告が出ないようなパラメタの設定(あるいはモデル式の修正)ができれば追記いたします。

 

4.結果の答え合わせ

ちょっとえげつない量の結果が出てくるのですが、とりあえずRhat<1.1は満たされているようです。
ちなみに、traceplotを見ると、やはりあまりきれいな結果にはなっていませんでした。
本当は修正が必要なのでしょうが、今回はそのまま進めます。

推定されたパラメタの答え合わせをします。

ちなみに、答えは以下の通りです。95%区間には一応入っているようです(標準偏差はすべて100で割らなければいけないことに注意してください)。
muZero:160
sigmaV:10(標準偏差です。シミュレーションの最後の段階でノイズを加えるときにsd=10と出てきています)
sigmaW:3(標準偏差です)
coefMedicineTrendZero:0.005
coefMedicineZero:-25
sigmaWcoef:0.5(標準偏差です)
sigmaWcoefTrend:0.03(標準偏差です)
sigmaVcoef:2(標準偏差です)

とりあえず、答えとそこまでかい離していなかったので、そのまま続けます。

5.推定結果の図示

時間的に変化する回帰係数を図示します。

まずは、係数の95%区間をとります。

後は、ggplot2を使って図示します。
実際の回帰係数(答え)を点線で表示するようにしました。

時変係数モデル_6_回帰係数プロット
(クリックすると大きくなります)

時間がたつにつれて、血圧が下がりにくい(回帰係数が0に近づく)ことがわかります。
答えも、一応95%信用区間には含まれているようです。

今回は若干荒いところもありましたが、複雑なモデルを推定する雰囲気が伝わればと思います。

パラメタの収束が悪かった件に関しては、パラメタの事前分布を変えるなど、いろいろのやり方が考えられます。
今後、追記、修正していく予定です。

 

参考文献

基礎からのベイズ統計学:
ハミルトニアンモンテカルロ法による実践的入門

 
何度も紹介していますが、ベイズ推定について書かれた良書です。
巻末の付録にStanの使い方も書かれています。
 

岩波データサイエンス Vol.1特集「ベイズ推論とMCMCのフリーソフト」


こちらも何度も紹介している本です。
ソフトの使い方に関しては、こちらのほうが詳しいです。
状態空間モデルは、季節調整モデルなどの解説があります。
 



スポンサードリンク