⏱️

一般化ベイズ法でハザード比が時間変化するCox回帰を実装する

2025/01/18に公開

この記事ではタイトルをStanで実装し、例題データであるBrem et al. (1995)に対して適用します。

一般化ベイズ法とCox回帰

菅澤先生の記事を参考にしました。ただし、コードの書き方は変えて、高速化を複数試しました。
https://qiita.com/ssugasawa/items/060dc8af0836ee209898
https://qiita.com/ssugasawa/items/e0f1ca1e19d23c643679

データ

データは以下の書籍のサポートページから入手できました

カプランマイヤー曲線は次の通りです。

ハザード比が時間変化しない通常のCox回帰

まずはハザード比が時間変化しない通常のCox回帰を次の2通りの方法で実装しました。

  • モデルa:log_sum_exp()を使って実装
  • モデルb:累積和を使って実装

Stanコード:モデル1a

https://github.com/MatsuuraKentaro/statistical_model_scraps/blob/main/cox_regression_with_time-varying_HR/stan/model1a.stan

  • 16~18行目:log_sum_exp()を使うと、ユニットごとの部分尤度を書く部分がこのように書けます。exp()を使わない分だけ安定と思われます。
  • 17行目:打ち切りのときはEventが0なので、targetへの寄与もゼロになります。
  • 22行目:リッジ回帰に相当します。

Stanコード:モデル1b

https://github.com/MatsuuraKentaro/statistical_model_scraps/blob/main/cox_regression_with_time-varying_HR/stan/model1b.stan

  • 14行目:cumulative_sum()reverse()をあわせて使うことで、ユニットごとの部分尤度を書く部分がベクトル化できます。その分高速になると思われます。ただしexp()は使うのでどうなるか。なお、次のRコードでデータを生存時間の逆順(長い順)に並びかえて渡せば、2つのreverse()を消すこともできます。
  • 19行目:リッジ回帰に相当します。

Rコード

https://github.com/MatsuuraKentaro/statistical_model_scraps/blob/main/cox_regression_with_time-varying_HR/run-model1.R

5行目:部分尤度の書き方から、データは生存時間の昇順(短い順)に並び替えておく必要があります。

結果

推定結果はほぼ完全に一致しました。計算時間はモデル1bがモデル1aに比べて約25倍高速でした。exp()とかは些細な問題で、ベクトル化するのが重要ということですね。以降ではモデル1bをベースにして時間変化版を実装します。

ハザード比が時間変化するCox回帰

生存時間が週数の実数であることも考慮し、ここではガウス過程を使って実装してみました。ガウス過程の部分は次の書籍の12章に従っています。
https://amzn.to/3E0xlaS

Stanコード:モデル2

https://github.com/MatsuuraKentaro/statistical_model_scraps/blob/main/cox_regression_with_time-varying_HR/stan/model2.stan

  • 62行目:リッジ回帰に相当します。

Rコード

https://github.com/MatsuuraKentaro/statistical_model_scraps/blob/main/cox_regression_with_time-varying_HR/run-model2.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