ねこすたっと

ねこの気持ちと統計について悩む筆者の備忘録的ページ。

生存時間データ分類(3):Competing type(tidycmprskパッケージ)[R]

生存時間データを、イベントの種類と1人あたり許容される回数をもとに4つに分けることは述べた。

necostat.hatenablog.jp

ここでは、「複数種類あるイベントのうち、いずれかを1回だけ経験しうる状況(競合型生存時間データ, competing type survival data)」について掘り下げたい。

例えば、

  • 担癌患者が、腫瘍が原因で死亡するか、心血管疾患が原因で死亡するか
  • 入院患者が、手術で死亡するか、院内感染症で死亡するか

などはこの状況にあたる。 このように「一方が他方の発生を阻害する関係にあるイベント」を競合リスク(competing risk)あるいは競合イベント(competing event)と呼ぶ。

図1:競合型生存時間データ

競合リスクを含んだ生存時間解析

競合リスクへの対処方法は

  1. 複合アウトカム(composite outcome):
    複数種のイベントを1つに統合し、いずれかのイベントを最初に発生した時点までの時間を解析に使う
  2. 原因別ハザードモデル(cause-specific hazard model):
    興味あるイベント以外のイベント発生については打ち切りとしてKaplan-Meier曲線を描く
  3. 部分分布ハザードモデル(subdistribution model):
    それぞれのイベントについて累積発生曲線(cumulative incidence function, CIF)を描く
  4. 多状態モデル(multi-state model):
    エントリーを含めたそれぞれの状態間の移行確率を考えて、ある時点で各状態に存在する確率を推定する。Aalen-Johansen推定法とも。

がある(と思う)。

1の代表は、心筋梗塞や脳卒中を複合アウトカムにしたmajor adverse cardiac event(MACE)だろう。1つにまとめてしまえば通常の生存時間解析と同じになる。 4は別の記事で考えてみる予定。なのでここでは2,3のみ扱う。

原因別ハザードモデル(cause-specific hazard model)

「興味あるイベント以外のイベントが発生した場合は打ち切りとして扱う」というシンプルな方法。 こうすれば、k番目*1のイベントに対するハザード関数  \lambda^{cs}_k(t) は「ある時点までどのイベントも経験しなかった人におけるk番目のイベントを瞬間発生率(rate)」と定義され、イベントの種類ごとに通常の方法でハザード関数が推定できる。

ただ問題になるのが、競合リスクが存在することにより、「打ち切り症例は打ち切りが発生した時点における生存集団を代表していなければならない」という独立性の仮定(independent censoring)を満たしていると考えられなくなること。こうなると、ハザード関数と累積発生関数(cumulative incidence function, CIF)は一対一対応でなくなる*2

もともとハザード(rate)は時間的局所において瞬間的な発生(risk)と直感的に結びつけて解釈されているので、CIFを適切に推定することができない方法は都合が悪い(実際、全てのイベントのCIFを足し合わせると1を超えてしまう)。累積発生数を適切にモデル化できていないので、予測を目的とした解析には適していない。

部分分布ハザードモデル(subdistribution hazard model)

Fine-Grayモデルとも呼ばれるこのモデルでは、k番目のイベントに対するハザード関数  \lambda^{sd}_k(t) は「ある時点までk番目のイベントを経験しなかった人におけるk番目のイベントを瞬間発生率(rate)」と定義される(競合リスクの存在のため、観察時間を無限大にしても累積発生関数(CIF)が1に到達しないため、厳密には分布関数とならないことが「部分分布(subdistribution)」という名称の由来)。

このモデルでは、イベントを全く経験していない人だけではなくて、その時点以前に競合イベントを発生したことがある人もat risk集団に含まれる。これは例えば、腫瘍関連死亡のハザードを考えるときに、心血管疾患が原因で死亡した人も対象集団に含まれていることを意味する。このため、因果推論を目的とした解析には適していない。

tidycmprskパッケージを使ってFine-Grayモデルに当てはめる

原因別ハザードモデルは、単に競合イベントを打ち切りとして扱うだけなので、通常の生存時間解析で使うcoxph( )を使えば良いとして、問題はFine-Grayモデルの方。

cmprskパッケージを使う方法と、timeregパッケージを使う方法がある。後者の方が解説は分かりやすかったが、もう少し一般的な累積発生割合をモデル化する方法と合わせて別の記事で。ここではcmprskパッケージをtidyverseと相性よく使えるtidycmprskパッケージを使う。

パッケージとデータセットの準備

tidycmprskパッケージに含まれるtrialというシミュレーションデータを使う*3

> library(tidycmprsk)
> head(trial)
# A tibble: 6 × 9
  trt      age marker stage grade response death death_cr           ttdeath
  <chr>  <dbl>  <dbl> <fct> <fct>    <int> <int> <fct>                <dbl>
1 Drug A    23  0.16  T1    II           0     0 censor                24  
2 Drug B     9  1.11  T2    I            1     0 censor                24  
3 Drug A    31  0.277 T1    II           0     0 censor                24  
4 Drug A    NA  2.07  T3    III          1     1 death other causes    17.6
5 Drug A    51  2.77  T4    III          1     1 death other causes    16.4
6 Drug B    39  0.613 T4    I            0     1 death from cancer     15.6
  • trt:化学療法
  • age:年齢
  • marker:腫瘍マーカー
  • stage:Tステージ
  • grade:グレード
  • response:腫瘍の反応性
  • death:死亡
  • death_cr:腫瘍関連死亡、腫瘍以外による死亡、打ち切り
  • ttdeath:死亡または打ち切りまでの時間

crr( )でモデルを当てはめる

書き方はcoxph( )と同じ感じ。

  • death_cr:イベントを示すfactor型の変数
  • failcode:上記の変数のカテゴリーの中で興味あるイベントを表すものを指定(指定しなかったら1番目のイベントが採用)
  • cencode:上記の変数のカテゴリーの中で打ち切りを表すものを指定
fit_crr <- crr(Surv(ttdeath, death_cr) ~ age + trt, 
               failcode = "death from cancer",
               cencode = "censor",
               data = trial_clean)

そのまま出力すれば表が出るが、指数変換はされていない。

> fit_crr

── crr() ───────────────────────────────────────────────────────────────────────
• Call Surv(ttdeath, death_cr) ~ age + trt
• Failure type of interest "death from cancer"

Variable    Coef    SE      HR     95% CI       p-value    
age         0.006   0.010   1.01   0.99, 1.03   0.56       
trtDrug B   0.417   0.279   1.52   0.88, 2.62   0.13    

せっかくなので、綺麗に要約した表を作ってくれるgtsummaryパッケージのtbl_regression( )と組み合わせてみる。gtsummaryパッケージの話はまた今度(まだまとめてないだけ)。

fit_crr %>% gtsummary::tbl_regression(exponentiate=TRUE)

表1:Fine-Grayモデルの結果をtbl_regression( )で綺麗に表示した

累積発生関数のグラフを描く

累積発生関数のグラフはcuminc( )とautoplot( )を組み合わせて。

cuminc(Surv(ttdeath, death_cr) ~ trt, 
       failcode="death from cancer",
       data = trial) %>%
  autoplot(conf.int = TRUE)

図2:cuminc( )とautoplot( )で累積発生関数のグラフを描いた

おわりに

  • いつも以上にコピペで使っただけ感が強いが、Fine-Grayモデルについてはtimeregパッケージをまとめるときにもう少し掘り下げたい。

参考資料

  • Therneau先生のRを使った生存時間解析のvignette: Multi-state models and competing risks

  • 競合リスクを考慮しないとどんな問題があるか、どのような方法で対処するのか、などが分かりやすく説明されているので、最初に読むのにおすすめ。

pubmed.ncbi.nlm.nih.gov

  • 競合リスクがあることでハザード関数(rate)と累積発生関数(risk)が一対一対応でなくなってしまうことを軸に説明してある。上を読んだ後にもう一回読んだら少し理解が進んだ(多分)

pubmed.ncbi.nlm.nih.gov

  • cmprskパッケージと結果をキレイに表示するtidycmprskパッケージの使い方(英語だけど分かりやすい)

cran.r-project.org

  • Aalen-Johnson推定について日本語で説明された記事。非常に詳しい。生存時間データ分類(4)でまとめるつもりだったけど、このページ読めばいいのでは?と思った。

tarohmaru.web.fc2.com

*1:発生した順序の話ではなく、K種類あるイベント中のk番目の種類ということ

*2:通常のKaplan-Meier推定法は独立性の仮定の下、両者が相互に変換可能

*3:cmprskパッケージには入ってないみたい