⏱️
一般化ベイズ法でハザード比が時間変化するCox回帰を実装する
この記事ではタイトルをStanで実装し、例題データであるBrem et al. (1995)に対して適用します。
一般化ベイズ法とCox回帰
菅澤先生の記事を参考にしました。ただし、コードの書き方は変えて、高速化を複数試しました。
データ
データは以下の書籍のサポートページから入手できました
カプランマイヤー曲線は次の通りです。
ハザード比が時間変化しない通常のCox回帰
まずはハザード比が時間変化しない通常のCox回帰を次の2通りの方法で実装しました。
- モデルa:
log_sum_exp()
を使って実装 - モデルb:累積和を使って実装
Stanコード:モデル1a
- 16~18行目:
log_sum_exp()
を使うと、ユニットごとの部分尤度を書く部分がこのように書けます。exp()
を使わない分だけ安定と思われます。 - 17行目:打ち切りのときは
Event
が0なので、target
への寄与もゼロになります。 - 22行目:リッジ回帰に相当します。
Stanコード:モデル1b
- 14行目:
cumulative_sum()
とreverse()
をあわせて使うことで、ユニットごとの部分尤度を書く部分がベクトル化できます。その分高速になると思われます。ただしexp()
は使うのでどうなるか。なお、次のRコードでデータを生存時間の逆順(長い順)に並びかえて渡せば、2つのreverse()
を消すこともできます。 - 19行目:リッジ回帰に相当します。
Rコード
5行目:部分尤度の書き方から、データは生存時間の昇順(短い順)に並び替えておく必要があります。
結果
推定結果はほぼ完全に一致しました。計算時間はモデル1bがモデル1aに比べて約25倍高速でした。exp()
とかは些細な問題で、ベクトル化するのが重要ということですね。以降ではモデル1bをベースにして時間変化版を実装します。
ハザード比が時間変化するCox回帰
生存時間が週数の実数であることも考慮し、ここではガウス過程を使って実装してみました。ガウス過程の部分は次の書籍の12章に従っています。
Stanコード:モデル2
- 62行目:リッジ回帰に相当します。
Rコード
- 9行目:予測する時点の数です。
- 10行目:予測する時点です。13行目で0週~200週を
に標準化して渡すので予測範囲も[0,1] になっています。[0,1]
結果
ハザード比の対数の時間変化をプロットしたものが以下になります。黒の実線は事後中央値、濃い灰帯は50%信用区間、薄い灰帯は90%信用区間です。
Enjoy!
- Brem, H. et al. (1995). Placebo-controlled trial of safety and efficacy of intraoperative controlled delivery by biodegradable polymers of chemotherapy for recurrent gliomas. The Lancet, 345(8956), 1008-1012.
Discussion