crispy.data

目次

    1. 生存時間解析とは
      1. 生存関数
      2. ハザード関数
      3. 生存関数とハザード関数の関係
    2. NCCTG Lung Cancer Data
      1. データの概観と分析の目的設定
      2. 打ち切りの確認
      3. ドメイン知識の確認
      4. データの確認
    3. 生存曲線の分析
      1. 推定生存曲線
      2. KM 曲線を対象にした統計検定

NCCTG Lung Cancer Dataで実践する生存時間解析

2023-01-28
エモリー大学本で学んだ内容を復習するため，実際のデータを使って生存時間解析をしてみる．

生存時間解析とは

エモリー大学本では，
興味のある結果変数（outcome variable）を，あるイベントが起こるまでの時間としたデータを解析する一連の統計的手法
とされている． すぐに思い付く例だと，
  • 治療を始めてから病気が治るまでの時間
  • 事故に遭ってから死亡するまでの時間 などが挙げられるが，人の生死に関する領域にしか適用できない訳ではなく，医学・生物学の分野以外でも活用できる手法である．
例えば，
  • サブスクサービスにおけるサブスクキャンセルまでの時間
  • マッチングアプリを始めてからマッチするまでの時間
  • 家電を使い始めてから買い換えるまでの時間 などを生存時間解析の枠組みで分析することも可能．
生存時間解析の枠組みでは，分析対象としたいこれらのイベントとイベント発生までの時間を生存関数やハザード関数としてモデリングし, 分析することが多い．

生存関数

生存関数 S(t)S(t) は，特定の時間 tt よりも長くイベントが発生しない確率を示す関数．表現がややこしいが, イベントを死亡と定義すれば, つまりは時間 tt 時点で生存している確率を示す関数. 具体的には．
S(0)=1S(∞)=0S(t+Δt)≤S(t)\begin{align*} & S(0) = 1 \\ & S(\infty) = 0 \\ & S(t+\Delta t) \leq S(t) \\ \end{align*}
の性質を持つ単調減少な関数．理論的には, 時間 tt によって滑らかに減少するが, 実際には観測する時間が離散的になるため, 非連続に変化する.

ハザード関数

ハザード関数は，単位時間 Δt\Delta t に対する 時点 tt まで生存した条件の下で時点 tt 〜 tt にイベントが発生する確率の比率である．ハザード関数自体はハザード"確率"ではなく，単位時間を基準にした相対的なリスク指標のようなもの．
h(t)=lim⁡Δ→0P(t≤T<t+Δt∣T≥t)Δth(t) = \lim_{\Delta \to 0} \frac{P(t \leq T \lt t + \Delta t | T \geq t)}{\Delta t}
ただし，
  • TT : 確率変数とした生存時間
  • t : TT の特定の値

生存関数とハザード関数の関係

生存関数とハザード関数は各時点の生存率に着目するか生存リスクに着目するかの違いだけであり, 一方が求まれば他方も求まる。
S(t)=exp[−∫0th(u)du] S(t) = exp\left[ - \int_0^t h(u) du \right]
h(t)=−[dS(t)/dtS(t)] h(t) = - \left[ \frac{d S(t) / dt}{S(t)} \right]
因みに, 臨床分野では生存リスク自体に感心がある場合があるからか, 生存確率を扱う生存関数よりもハザード関数の推定やハザードを目的としたモデリング手法が発展している印象がある. (まだまだ入門したてのため, 誤認かもしれない.) 生存時間解析ではこれらの関数を推定することで，イベントが発生しやすい時期の特定や複数の群の中で特にイベントが早期に発生しやすい群の特定などを行うことができる．

NCCTG Lung Cancer Data

生存時間解析の概要を示したところで，今回分析したいデータを紹介する.

データの概観と分析の目的設定

本投稿では survivalパッケージに入っている NCCTG Lung Cancer Data を扱う． NCCTG Lung Cancer Dataは，進行肺がん患者の生存期間と、医師および患者自身によって測定された患者のパフォーマンスステータスの評価が記録されている。 早速データの中身を見てみる．
library(tidyverse)
library(survival)
lung %>% head()
＃insttimestatusagesexph.ecogph.karnopat.karnomeal.calwt.loss
1330627411901001175NA
23455268109090122515
331010156109090NA15
45210257119060115011
518832601010090NA0
61210221741150805130
下記 10 項目の情報が含まれており，患者一人につき 1 行のデータが与えられている形式らしい． inst: 施設コード time: 生存日数 status: 1 = 打ち切り, 2 = 死亡 age: 年齢 sex: 1 = 男性, 2 = 女性 ph.ecog: 医師が評価した ECOG パフォーマンススコア. 0 = 無症状, 1 = 症状はあるが歩行可能な状態, 2 = 半日未満ベッドにいる状態, 3 = 寝たきりではないが半日以上ベッドにいる状態, 4 = 寝たきりの状態 ph.karno: 医師が評価した Karnofsky パフォーマンススコア pat.karno: 患者が評価した Karnofsky パフォーマンススコア meal.cal: 食事で摂取したカロリー wt.loss: 過去 6 ヶ月での体重減少量（ポンド） 本データセットを使った生存時間解析の目的はいろいろ考えられそうだが，今回はph.karnoとpat.karnoのスコアによって患者の生存リスクに差異が生じているか，及びどちらがより生存リスクを推定する共変量として有効であるかについて分析することを目標にする．

打ち切りの確認

statusに記載している"打ち切り"は，特定期間のデータが観測できていないような状態を示す．ここでは打ち切りの代表的な類型として，右側打ち切りと左側打ち切りを図で説明する．
図中 A は打ち切りが発生していないケースである．開始時点から at-risk な状態が確認されており，その後観測期間中にイベントが発生している． 図中 B は途中から実験対象が追跡できなくなってしまうなどの右側打ち切りのケースである．観測期間中にイベントが発生していない場合は全て右側打ち切りになるため，打ち切りの中でも右側打ち切りが最も多いと思われる． 図中 C はイベントのみが観測されており，それ以前の観測したデータが存在しないケースである左側打ち切りを表している．例えば感染症への感染から症状発生までを扱う場合において，0 日目から 10 日目の間に感染したことは明らかだが，10 日時点の症状発生しか観測できていないような場合がこれに該当する． 打ち切りの概念が存在することで，途中で観測できない期間が生じているような個体のデータも含めて分析対象にできるというメリットがある． 例えば上図の例だと，10 日の時点でat-riskな状態であった B と C にイベントが発生したことが確認でき，10 日を境に生存率が激減（100% → 33%）する現象をより精緻に確認できる．（打ち切りのデータを除くと，100% → 0%の生存率の変化しか扱えず，実態に合わない．） ちなみに本データセットでのstatusの説明に記載されている"打ち切り"は，time 以降の観測データが存在しないことを示しているため，具体的には右側打ち切りという状態を表現していることになる．

ドメイン知識の確認

ECOG パフォーマンススコアと Karnofsky パフォーマンススコアについて知らなかったので調べてみた．どちらも症状の深刻さを表す点数らしい． なおデータを見て分かるように, ECOG スコアは医師が評価したもののみが存在するが, Karnofsky スコアは医師が評価したものと患者が自己評価したものの 2 つが存在する. ECOG パフォーマンススコアの定義（http://www.jcog.jp/doctor/tool/C_150_0050.pdf より）
スコア定義
0全く問題なく活動できる。 発病前と同じ日常生活が制限なく行える。
1肉体的に激しい活動は制限されるが、歩行可能で、軽作業や座っての作業は 行うことができる。 例：軽い家事、事務作業
2歩行可能で自分の身の回りのことはすべて可能だが作業はできない。 日中の 50% 以上はベッド外で過ごす。
3限られた自分の身の回りのことしかできない。日中の 50% 以上をベッドか 椅子で過ごす。
4全く動けない。 自分の身の回りのことは全くできない。 完全にベッドか椅子で過ごす。
Karnofsky パフォーマンススコアの定義（http://www.npcrc.org/files/news/karnofsky_performance_scale.pdf より）
概要スコア定義
日常生活や仕事をすることができる．特段ケアが必要ない状態．100特に不調がなく，通常の体調．病気である兆候も存在しない．
90日常生活を送れる．病気の兆候や症状は軽微である．
80日常生活を送るために努力を要する．病気の兆候や症状が存在する．
仕事はできないが自宅で生活を送ることができ，ほとんどの身の回りのことを自身で行うことができる．必要な介助の程度は場合によって変わる．70自活しているが，日常生活や活動的な仕事はできない．
60時折介助が必要だが，身の回りのことはほとんど自分でできる．
50かなりの介助が必要であり，頻繁に医療行為が必要になる．
身の回りのことを自身で行えない．施設や病院の介助に相当するものが必要．病状が急速に悪化することもある．40障害．特別なサポートや介助を必要とする。
30重度障害．死の危険はないが、入院が必要な状態。
20重症．入院が必要．積極的な支持療法が必要。
10瀕死．致命的な状況が急速に進行している。
0死亡．
Karnofsky パフォーマンススコアは既に cardinality が低いデータだが，3 水準にまとめ上げるような使い方ができそう．

データの確認

lung %>%
  mutate(status = ifelse(status == 1, "censored", "dead")) %>%
  group_by(status) %>%
  summarise(n = n()) %>%
  ggplot(aes(x = status, y = n)) +
  geom_bar(stat = "identity") +
  ylab(element_blank())
大体1:3ぐらいの比率でdeadのデータが多い． Karnofskyスコアの定義で確認したように，ドメイン知識を基に3水準に分けれそうなので，040, 5070, 80~100に分けて度数を確認してみる．
lung <-
  lung %>%
  mutate(
  binned_ph.karno = case_when(
    ph.karno >= 0 & ph.karno <= 40 ~ "0~40",
    ph.karno >= 50 & ph.karno <= 70 ~ "50~70",
    ph.karno >= 80 & ph.karno <= 100 ~ "80~100",
    TRUE ~ "missing"
  ),
  binned_pat.karno = case_when(
    pat.karno >= 0 & pat.karno <= 40 ~ "0~40",
    pat.karno >= 50 & pat.karno <= 70 ~ "50~70",
    pat.karno >= 80 & pat.karno <= 100 ~ "80~100",
    TRUE ~ "missing"
  )
)

lung %>%
  ggplot(aes(x = binned_ph.karno)) +
  geom_bar()

lung %>%
  ggplot(aes(x = binned_pat.karno)) +
  geom_bar()
どちらもごく一部のデータが欠損しているらしい． また，binned_ph.karno （医師が評価したKarnofslyパフォーマンススコア）に関しては，0~40をつけた例が無いことが確認できる．
lung %>%
  ggplot(aes(x = pat.karno, y = ph.karno)) +
  geom_count() +
  geom_smooth(method = "lm", formula = y ~ x)
  • pat.karno が 80 未満の場合は，医師の診察結果よりも自覚症状の方が重い傾向にある．
  • pat.karno が 80 以上の場合は，自覚症状よりも医師の診察結果の方が重い傾向にある．

生存曲線の分析

分析の目的を考え，karnoの値が欠損しているレコードは削除し，Karnofskyパフォーマンススコアのまとめ上げは，070，80100 の 2 水準にまとめ上げることにする．
lung <-
  lung %>%
  mutate(
    binned_ph.karno = case_when(
      ph.karno >= 0 & ph.karno <= 70 ~ "0~70",
      ph.karno >= 80 & ph.karno <= 100 ~ "80~100",
      TRUE ~ "missing"
    ),
    binned_pat.karno = case_when(
      pat.karno >= 0 & pat.karno <= 70 ~ "0~70",
      pat.karno >= 80 & pat.karno <= 100 ~ "80~100",
      TRUE ~ "missing"
    )
  ) %>%
  filter(
    binned_ph.karno != "missing",
    binned_pat.karno != "missing"
  )

推定生存曲線

カプランマイヤー曲線（以下，KM 曲線）は，X 軸が経過時間，Y 軸が推定生存率を表すような可視化手法． ここで推定生存率は，
S~(t)=S~(t−1)×nt−ftnt\tilde{S}(t) = \tilde{S}(t-1) \times \frac{n_t - f_t}{n_t}
式から明らかなように, 再発や観測途中から at-risk なサンプルが増えるような例外的なケースを除いて，原則として単調現象な階段状のグラフになる． とりあえず，2水準にまとめ上げたph.karno と pat.karno のKM曲線を出力してみる．
library(ggsurvfit)

survfit(
  Surv(time, status)~binned_ph.karno,
  data=lung
) %>%
  ggsurvfit() +
  xlab("days")

survfit(
  Surv(time, status)~binned_pat.karno,
  data=lung
) %>%
  ggsurvfit() +
  xlab("days")
患者が評価したか医師が評価したかにかかわらず，Karnofsky が高い群（80 ~ 100）の方が生存率が高い傾向にあることが確認できる．症状が軽い方が生存率が高いということであり，当然の結果になっている． 特筆すべきは，binned_ph.karnoのKM曲線について，800 日以降ぐらいの傾向を確認すると生存率が逆転しているという点．医師が評価した Karnofsky スコアよりも患者が評価した Karnofsky スコアの方が信頼できる(≒ 病状の深刻さを反映している)のではないかという仮説が考えられる． しかし，そもそも Karnofsly スコアが高い群と低い群で生存率が異なる傾向が確認できているものの，KM 曲線を眺めているだけでは偶然生じた差異であることを否定できない．KM 曲線の差が統計的に有意であるかを確認してみる．

KM 曲線を対象にした統計検定

KM 曲線に有意差が存在するかを確認する手法の内 5 つを紹介する．他の手法も含めて内容を詳しく知りたい方は，こちらのペーパーを見るといい．