Top 5 This Week

関連記事

【医療統計版】1. 生存時間解析入門:打ち切りがあるデータで何を推定するのか

- 本サイト運営者のサービスの紹介 -

多施設の観察研究のデータを渡されて、共著者からこう言われたことはないでしょうか。「追跡できなくなった12例は除外して集計してほしい」。転居や来院中断で最後まで追えなかった患者を落として、残った症例だけで平均生存期間を出す。手順としては単純で、集計表もきれいに揃います。

ただ、この操作は生存期間を系統的に短く見せます。しかも標本を増やしても改善しません。

生存時間解析は、この「最後まで追えなかった人のデータをどう使うか」という一点のために組み立てられた枠組みです。ここでは打ち切りを捨てたときに推定値がどれだけ動くかを実際に計算し、そのうえで生存関数とハザード関数が何を分担しているのか、分布の仮定をどこまで置くべきかを見ていきます。

生存時間解析が答えている問い

生存時間解析が対象にするのは、ある事象が起きるまでの時間です。「何人がイベントを経験したか」ではなく「いつ起きるか」を推定の対象に置きます。この違いは、そのまま扱えるデータの形の違いになります。

個体iについて記録するのは、観察開始からの時間$t_i \ge 0$と、イベントが実際に起きたかどうかを表す指示変数$\delta_i \in \{0,1\}$の組です。$\delta_i = 1$ならその時点でイベントが起きたこと、$\delta_i = 0$なら観察を終えた時点でまだ起きていなかったことを表します。1人あたり数値が2つ必要になる、というのがこのデータ構造の要点です。

ロジスティック回帰との違いは、この2つ目の数値を使うかどうかにあります。ロジスティック回帰は「5年以内に再発したか」を0と1に潰してしまうので、4年11か月で再発した患者と、4か月で再発した患者が同じ扱いになります。3年で追跡が切れた患者に至っては、0と1のどちらにも置けません。

比較軸 線形回帰 ロジスティック回帰 生存時間解析
従属変数 連続量 二値 非負の時間と発生指示変数の組
打ち切りの扱い 枠組みを持たない 枠組みを持たない 尤度への寄与を分けて扱う
答える問い 期待値はいくつか 発生確率はいくつか いつ起きるか
代表的手法 最小二乗法 最尤法 Kaplan-Meier法、Cox回帰

イベントの定義は1つに絞ります。「死亡または再発」を主要評価項目にすると、どちらが起きたのかで解釈が変わるため、複合エンドポイントとして扱うか競合リスクとして扱うかを先に決める必要があります。定義が二重になったまま解析に入ると、あとから直せません。

打ち切りを捨てると何が起きるか

右打ち切りは、観察を終えた時点でまだイベントが起きていない状態を指します。試験終了まで生存していた患者、途中で来院しなくなった患者、他院に移った患者。臨床研究のデータでは例外ではなく常態です。

無作為化からの追跡期間を6名分並べた図。3名はイベント発生、3名は打ち切りで、うち1名は試験終了時点で打ち切られている。
図1 6名分の追跡。P4は試験終了時点の45か月で打ち切られており、これは研究計画によって生じる打ち切りです。一方P2とP6は追跡が途切れた時点で打ち切られています。同じ記号でも発生の理由が違う点が、あとで無情報打ち切りを考えるときに効いてきます。

打ち切りの観測は情報を持っていないわけではありません。P2は32か月時点まではイベントを起こさなかった、という事実が確定しています。これは「32か月までは生存した」という制約であり、生存関数の推定に使えます。

では、この情報を捨てるとどうなるか。真の生存時間の中央値が24か月と分かっている架空データを作り、全例を使ったKaplan-Meier推定と、打ち切り例を除いた推定を比べてみます。以下は説明のための架空データです。

library(survival)

set.seed(42)
n <- 600
true_time <- rexp(n, rate = log(2) / 24)   # 真の生存時間、中央値は24か月
cens_time <- runif(n, min = 0, max = 60)   # 追跡が切れる時点
obs_time  <- pmin(true_time, cens_time)
status    <- as.numeric(true_time <= cens_time)

km   <- survfit(Surv(obs_time, status) ~ 1)
drop <- survfit(Surv(obs_time[status == 1], status[status == 1]) ~ 1)

summary(km,   times = 24)
summary(drop, times = 24)
#>  time n.risk n.event survival std.err lower 95% CI upper 95% CI
#>    24    200     218    0.553  0.0233         0.51        0.601
#>
#>  time n.risk n.event survival std.err lower 95% CI upper 95% CI
#>    24     65     218     0.23   0.025        0.186        0.284

24か月時点の真の生存率は0.500です。全例を使った推定は0.553、95%信頼区間は0.51から0.60でした。打ち切り例を除いた推定は0.23です。

ここは少し注意が必要です。全例を使った0.553も真値からずれており、信頼区間の下限がぎりぎり0.500を上回っています。標本を1つ取った以上、このくらいのずれは起こります。区別したいのは、ずれの性質のほうです。0.553は別の乱数を引けば反対側にもずれます。0.23は何回引いても下にずれます。

真の生存関数と、600例すべてを使ったKaplan-Meier推定、打ち切り例を除外した推定の3本を重ねた図。打ち切り例を除いた曲線だけが大きく下に外れている。
図2 破線が真の生存関数、実線2本が推定値です。全例を使った曲線は破線に沿って動きますが、打ち切り例を除いた曲線は追跡が進むほど下へ離れていきます。除外した患者は「まだイベントが起きていない人」なので、残った標本はイベントが早く起きた人に偏ります。
plot(km, conf.int = FALSE, xlab = "Months", ylab = "S(t)")
lines(drop, conf.int = FALSE)
curve(exp(-log(2) * x / 24), from = 0, to = 60, lty = 2, add = TRUE)

この偏りが尤度のどこから来るのかは、式にすると1行で見えます。イベントが観測された個体は密度$f(t_i)$を、打ち切られた個体は生存関数$S(t_i)$を尤度に寄与させます。

$$L = \prod_{i=1}^{n} f(t_i)^{\delta_i} \, S(t_i)^{1-\delta_i} \tag{1}$$

式(1)は、イベントが起きた人からは「ちょうどその時刻に起きた」という情報を、打ち切られた人からは「少なくともその時刻までは起きなかった」という情報を、それぞれ別の形で取り出していることを表しています。打ち切り例を除外するというのは、右側の因子をまるごと捨てるのと同じです。

ここまでは、打ち切りがイベントの起きやすさと無関係に生じるという前提で話してきました。これを無情報打ち切りと呼びます。有害事象が重くなった患者ほど脱落しやすい、といった状況ではこの前提が崩れ、Kaplan-Meier推定量も一致性を失います。脱落理由が判明している場合は理由別の集計を論文に載せることを勧めます。理由が全例「不明」と書かれた表は、査読者が最初に指摘する箇所です。

生存関数とハザード関数の役割分担

生存時間解析には、同じ分布を別の角度から見る2つの関数が出てきます。どちらか一方で足りるのに2つ用意されているように見えて、実際には答えている質問が違います。

生存関数$S(t)$は、時刻$t$までイベントが起きない確率です。

$$S(t) = \Pr(T > t) \tag{2}$$

ここで$T$はイベント発生時刻を表す確率変数です。式(2)は「今から$t$まで持ちこたえる確率」を表しており、$S(0)=1$から始まって単調に減っていきます。患者に説明するときに使えるのは、ほぼこちらです。

ハザード関数$h(t)$のほうは、時刻$t$まで生存した人が、その直後の短い時間にイベントを起こす瞬間的なリスクを表します。

$$h(t) = \lim_{\Delta t \to 0} \frac{\Pr(t \le T < t + \Delta t \mid T \ge t)}{\Delta t} \tag{3}$$

式(3)の条件付き確率の部分が重要で、分母の条件は「時刻$t$の時点でまだイベントを起こしていない」ことです。つまりハザードは、その時点で残っている人だけを見た危険度になります。術後1か月の死亡ハザードが高いのは、術後1か月まで生きている人にとってのリスクの話であって、患者全体に占める割合ではありません。

この2つは互いに変換できます。ただ、モデルを組み立てるときはハザードの側から入るほうが自然です。共変量の効果を「ハザードを何倍にするか」で表せば、時間の関数である基準ハザードを特定しないまま係数を推定できるからです。Cox比例ハザードモデルが広く使われている理由はここにあります。関数どうしの変換式と導出は生存関数とハザード関数の回で扱います。

中央値がすべて24か月で一致する3本の生存曲線。ハザード一定の指数分布、ハザード増加のワイブル分布、ハザード減少のワイブル分布で、中央値の前後の形が大きく異なる。
図3 3本とも中央値は24か月で一致していますが、ハザードの形が違うため曲線の形が違います。中央値だけを比較すると、この3つは区別がつきません。生存期間中央値を1つ報告して終わりにすると何が落ちるのか、この図が答えになっています。

分布の仮定をどこまで置くか

解析手法は、生存時間の分布にどこまで踏み込むかで分かれます。仮定を置かないほど柔軟ですが、そのぶん推定は不安定になります。

分布形を仮定しないのがノンパラメトリック推定で、Kaplan-Meier推定量とNelson-Aalen推定量がこれにあたります。観察された時点の集合の上でしか値が決まらないため、追跡期間の外へ外挿することはできません。5年追跡の研究から10年生存率を読むことはできない、というのはここから来ています。

共変量の効果だけをパラメトリックに推定し、基準ハザードは形を決めないままにするのがセミパラメトリックなアプローチです。Cox回帰がその代表で、医学論文の多変量解析では事実上の標準になっています。

分布を完全に指定するのがパラメトリックモデルです。指数分布、ワイブル分布、対数正規分布などを当てはめます。指定が正しければ推定効率は最も高くなりますが、外れたときは推定値が真値に近づきません。

どれを使うかで迷ったとき、私ならまずKaplan-Meier曲線を描いて研究の主目的を確認します。群間で生存を比べたいだけならKaplan-Meier法とログランク検定で足り、共変量の調整が必要ならCox回帰に進みます。ここでいきなりワイブル分布を当てはめる進め方は勧めません。分布形の当てはまりを論文で議論する負担が増えるわりに、ハザード比という共通言語での報告から遠ざかるからです。

ただし、追跡期間を超えた時点の生存率を推定する、あるいは医療経済評価で生涯の期待値を出す、といった外挿が目的に含まれるなら話が変わります。この場合は最初からパラメトリックモデルを選び、分布形の選択根拠と感度分析を本文に書くことになります。手法選択は分布の性質ではなく、何を報告する必要があるかで決まります。パラメトリック生存モデルの回では、この当てはまりの評価を扱っています。

論文でどう報告するか

報告ガイドラインは、生存時間解析について具体的な記載を求めています。観察研究についてのSTROBE声明(von Elm et al., Ann Intern Med 2007;147:573-577)は、追跡期間の要約と各群の脱落数を項目として挙げています。無作為化比較試験のCONSORT声明も、追跡と除外の流れを図で示すことを求めています。

実際の投稿で指摘が入りやすいのは、追跡期間の書き方です。「追跡期間中央値36か月」とだけ書かれていても、それが観察時間の中央値なのか、逆Kaplan-Meier法で求めた値なのかで意味が変わります。イベントが多い研究では前者が短めに出るため、Schemper and Smith(Control Clin Trials 1996;17:343-346)は逆Kaplan-Meier法を勧めています。

もう1つ、生存期間中央値に到達していないときの書き方があります。追跡期間内に生存率が0.5を下回らなければ中央値は推定できません。ここを空欄にしたり、最終観察時点の値で埋めたりせず、「中央値未到達(NR)」と書いたうえで、24か月生存率のような特定時点の推定値を信頼区間つきで併記します。この形にしておくと、あとから追跡を延長した続報とも比較できます。

この枠組みが作られてきた道筋

打ち切りを扱う発想そのものは、19世紀の保険数理の生命表にさかのぼります。ただ、現在の臨床研究で使う道具が揃ったのは20世紀後半です。

1958年、KaplanとMeierがJournal of the American Statistical Associationに積極限推定量を発表しました。分布形を仮定せずに打ち切りデータから生存関数を推定できるこの方法は、そのまま今も使われています(Kaplan-Meier推定量)。

1972年のCoxによる比例ハザードモデルは、基準ハザードを特定せずに共変量の効果を推定する道を開きました。部分尤度という推定の枠組みが正当化されるまでには議論があり、その理論的整理は1970年代を通じて進みます。2010年代以降はRandom Survival Forestのような機械学習側の手法も加わりましたが、予測性能と引き換えに、効果量の解釈と仮定の透明性という点では従来手法と違う課題を抱えています。

次に何を確かめるか

自分のデータで最初に見るべきは、打ち切りの数と理由です。全体の何割が打ち切りで、そのうち試験終了によるものと追跡中断によるものがどれくらいか。この2つの比率が分かれば、無情報打ち切りの前提がどれだけ危ういかの見当がつきます。そのうえでKaplan-Meier曲線を描き、中央値だけでなく特定時点の生存率も見ておく。ここまでが、どの手法に進むかを決める前の作業になります。

打ち切りと切断の区別、そして左切断が入ったときの扱いは、打ち切りと切断の回で詳しく見ます。

ここまで読んで、自分の研究では打ち切りの割合が高すぎるのではないか、脱落理由をどう書けばよいのか、と気になった方もいるかもしれません。実際のデータでは、打ち切りの発生時期、脱落理由の分布、事前に定めた解析計画のどこまでを本文に書くかまで含めて判断することになります。Dr.データサイエンスでは、こうした解析方針と報告のご相談を承っています。

参考文献

Kaplan EL, Meier P. Nonparametric estimation from incomplete observations. Journal of the American Statistical Association. 1958;53(282):457-481.

Cox DR. Regression models and life-tables. Journal of the Royal Statistical Society: Series B. 1972;34(2):187-202.

Schemper M, Smith TL. A note on quantifying follow-up in studies of failure time. Controlled Clinical Trials. 1996;17(4):343-346.

von Elm E, Altman DG, Egger M, et al. The Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) statement. Annals of Internal Medicine. 2007;147(8):573-577.

Klein JP, Moeschberger ML. Survival Analysis: Techniques for Censored and Truncated Data. 2nd ed. Springer; 2003.

Popular Articles