[Stan] 二つの時系列データの間に「差」があるか判断するには
2016-02-19
詳しい経緯はこのまとめを参照してください。時間軸でぶった切って各時点で検定を使う手法は、百歩譲って「差があるかどうか」は判定できるかもしれないけど、「どれほど異なるのか」については何も言えない。「どの時刻から異なるか」についても言えるか分からない。そこでベイズ統計モデリングで判断しようと思います。ベイズ統計モデリングでは多くの事前知識を仮定としてモデルに組み込みますが、検定も多くの仮定を前提にしている点は同様と思います。データは雰囲気だけ似せて自作しました。野生型100個体、変異体10個体で1～24まで1時間ずつ測定して24時点としました。まとめを見ると144時間みたいですが24時間に簡略化します。データの構成は以下です。
| type | X1 | X2 | … | X23 | X24 |
|---|---|---|---|---|---|
| 0 | 0.071 | 0.555 | … | -0.236 | -0.597 |
| 0 | 0.445 | 0.483 | … | -0.149 | 0.231 |
| 0 | 0.225 | 0.764 | … | -0.116 | -0.213 |
| … | … | … | … | … | … |
| 1 | 0.28 | 0.379 | … | -0.255 | -0.289 |
| 1 | 0.756 | 0.432 | … | -0.592 | 0.162 |
| 1 | 0.493 | 0.142 | … | -0.437 | -0.038 |
平均と±標準偏差のエラーバーつきで図にすると以下になります。
●モデル1
はじめに考えたモデルは以下です。
モデル式にしますと以下になります。各時刻で野生型と変異体と同じ状態空間モデルに従うベースラインを持つとする。変異体はそれに加えて定数の項（または指数関数で落ちる項）を付け加えて、変化点検出込みでモデリング。推定後、定数の高さのベイズ信頼区間から判定、でよさそう。 https://t.co/Aq0WDCC6P3
— .(@berobero11) 2016, 2月 16
元データは綺麗なsinカーブでしたが、ここでは一般の時系列でも通用する状態空間モデルを使っています。ここで、
状態空間モデルを使った推定については岩波データサイエンスvol.1が充実しています。また、この変化点検出のモデルは、BUGS Book 11.7.1節を参考にしています。
Stanで実装すると以下になります。
・3～4行目：「0」がついている方が野生型、「1」がついている方が変異体です。
・22行目：Stanではこのようにif_else文を使うことができます。
・29行目：for (t in 1:N_time) Y0[n,t] ~ normal(mu[t], s_Y); と同じ意味です。Stanのvector型を使うとこのように短縮形で書くことができ、関数を呼ぶ回数が減るので実行速度が速くなります。
・31行目：29行目と同様です。
推定結果は以下になります。diffの95%ベイズ信頼区間が0を含んでおらず、野生型と変異体の差はあると言えるでしょう。
| mean | se_mean | sd | 2.5% | 25% | 50% | 75% | 97.5% | n_eff | Rhat | |
|---|---|---|---|---|---|---|---|---|---|---|
| mu[1] | 0.23 | 0.00 | 0.03 | 0.17 | 0.21 | 0.23 | 0.25 | 0.29 | 887 | 1.01 |
| mu[2] | 0.49 | 0.00 | 0.03 | 0.44 | 0.48 | 0.49 | 0.51 | 0.55 | 843 | 1.01 |
| … | … | … | … | … | … | … | … | … | … | … |
| mu[23] | -0.27 | 0.00 | 0.03 | -0.33 | -0.29 | -0.27 | -0.25 | -0.22 | 1107 | 1.00 |
| mu[24] | 0.04 | 0.00 | 0.03 | -0.02 | 0.02 | 0.04 | 0.06 | 0.10 | 1151 | 1.00 |
| s_mu | 0.20 | 0.00 | 0.03 | 0.15 | 0.18 | 0.20 | 0.22 | 0.28 | 401 | 1.00 |
| s_Y | 0.31 | 0.00 | 0.00 | 0.30 | 0.31 | 0.31 | 0.31 | 0.32 | 1254 | 1.00 |
| diff | 0.32 | 0.00 | 0.05 | 0.23 | 0.29 | 0.32 | 0.36 | 0.42 | 174 | 1.01 |
| cp1 | 6.45 | 0.01 | 0.35 | 5.56 | 6.21 | 6.47 | 6.73 | 6.97 | 568 | 1.00 |
| period | 5.69 | 0.04 | 0.74 | 4.30 | 5.16 | 5.71 | 6.20 | 7.15 | 323 | 1.01 |
| y1_mean[1] | 0.23 | 0.00 | 0.03 | 0.17 | 0.21 | 0.23 | 0.25 | 0.29 | 887 | 1.01 |
| … | … | … | … | … | … | … | … | … | … | … |
| y1_mean[24] | 0.04 | 0.00 | 0.03 | -0.02 | 0.02 | 0.04 | 0.06 | 0.10 | 1151 | 1.00 |
| lp__ | 1798.55 | 0.13 | 3.78 | 1790.09 | 1796.25 | 1798.90 | 1801.24 | 1804.94 | 797 | 1.00 |
共通のベースライン（mu；左図）と、変異体と野生型の差（右図）を描画すると以下になります。
ここで、黒点とつなぐ線はMCMCサンプルの中央値、濃い灰帯は50%ベイズ信頼区間、薄い灰帯は95%ベイズ信頼区間です。
このモデルでは差を「一定期間の定数」と仮定していますが、定数じゃないかもしれません。また、変化点の外側でも差がある可能性を残す方が自然に思えます。そして、このようなif文を使ったMCMCが非効率的で初期値の設定もシビアな上に計算に時間を要します（1chainあたり5分ぐらいですが）。そこで、モデル2を考えました。
●モデル2
モデル式は以下です。
システムモデルの変化がなめらかであるとして、2階差分に変更しました（ちなみにモデル1で2階差分に変更すると、さらに推定に時間がかかります）。そして、野生型と変異体の差の部分も時間変化するようにし、変化点検出を考慮してシステムモデルにコーシー分布を使いました。
・13,15,24行目など：コーシー分布を使った状態変化はそのまま実装すると収束しません。一様分布に従うパラメータから変数変換で作るのがポイントです。詳しくは以前の記事を参照してください。
・25行目：for (t in 1:N_time) y1_mean[t] <- mu[t] + mu[di];をvector型の活用によって短縮しています。
・30行目：2階差分にしています。
こちらは計算時間は1chainあたり約15秒でした。推定結果は以下の通り。
共通のベースライン（mu；左図）と、変異体と野生型の差（di；右図）を描画すると以下になります。凡例は先ほどの図と同じです。
diの95%ベイズ信頼区間が0を含んでいない期間が差がある期間と言えるでしょう。さらに、どこから差がありそうなのか、どれほど差がありそうなのかも確率付きで述べることができます。
最後に、参考までに実行するRコードは以下になります。
・16行目：シビアな初期値設定に注目してほしいです。
・17行目：モデル1ではMCMCの自己相関が強く、iterを4500にしています。
データを作成したRコードは以下になります。
- 関連記事