投稿した原稿のResultsに「5年生存率は62%（95% CI 41-78）でした」と書いたとします。査読者から返ってきたのは、こういうコメントです。「Please add the number at risk below the survival curves.」
リスクテーブルを作ってみると、5年時点で追跡中だったのは7人でした。62%という数字は、この7人の上に載っています。
Kaplan-Meier推定量は打ち切りを含むデータから生存関数を推定してくれますが、時点によって推定の確からしさが大きく違います。ここでは推定量が積の形になる理由を確認したうえで、信頼区間の作り方で右端の見え方がどう変わるか、そして曲線のどこまでを結論の根拠にしてよいかを見ていきます。
なぜ掛け算になるのか
イベントが観測された時刻を$t_{(1)} < t_{(2)} < \cdots$と並べ、時刻$t_{(j)}$の直前まで追跡できていた人数を$n_j$、そのときのイベント数を$d_j$とします。この時刻を「越えられる」条件付き確率は$1 – d_j/n_j$です。
時刻$t$までの生存確率は、それまでの各時点を順に越えてきた確率の積になります。
$$\hat{S}(t) = \prod_{t_{(j)} \le t} \left(1 – \frac{d_j}{n_j}\right) \tag{1}$$
式(1)は、生存曲線がイベントのあった時点でだけ下に跳び、それ以外では平らなままになることを表しています。打ち切りの時点では新しい項が入らないので値は動きませんが、次のイベント時点での$n_j$が小さくなります。打ち切りは曲線を下げるのではなく、その後の1回ごとの下げ幅を大きくする、という効き方をします。
10例の架空データで確かめます。観測時間は3, 4, 5, 5, 8, 8, 9, 11, 12, 14か月で、このうちイベントは3, 5, 8, 8, 11か月の5件、残りは打ち切りです。以下は説明のための架空データです。
library(survival)
months <- c(3, 4, 5, 5, 8, 8, 9, 11, 12, 14)
event <- c(1, 0, 1, 0, 1, 1, 0, 1, 0, 0)
fit <- survfit(Surv(months, event) ~ 1)
summary(fit)#> time n.risk n.event survival std.err lower 95% CI upper 95% CI
#> 3 10 1 0.900 0.0949 0.732 1.000
#> 5 8 1 0.787 0.1340 0.564 1.000
#> 8 6 2 0.525 0.1759 0.272 1.000
#> 11 3 1 0.350 0.1849 0.124 0.985
4番目のイベントでは、たった1件で曲線が0.525から0.350へ、17.5ポイント落ちています。同じ1件でも最初は10ポイントでした。曲線の右側で段差が大きくなるのは、予後が急に悪くなったからではなく、分母が小さいからです。
信頼区間の作り方で右端の見え方が変わる
上の出力で気になる列があります。95%信頼区間の上限が、3つの時点で1.000ちょうどになっています。生存率が1という上限は、統計的な情報というより計算の打ち切りです。
分散はGreenwoodの公式で推定します。
$$\widehat{\mathrm{Var}}[\hat{S}(t)] = \hat{S}(t)^2 \sum_{t_{(j)} \le t} \frac{d_j}{n_j (n_j – d_j)} \tag{2}$$
式(2)は、各イベント時点での不確実性が積み上がっていく形になっています。$n_j$が小さいほど1項あたりの寄与が大きく、右へ行くほど分散は増えます。
問題はここから先で、この分散を使って区間をどう作るかに3通りあります。推定値に標準誤差の1.96倍を足し引きする線形の方法、対数変換してから戻す方法、そして$\log(-\log \hat{S}(t))$を変換して戻す方法です。Rの既定は2番目です。
for (type in c("plain", "log", "log-log")) {
f <- survfit(Surv(months, event) ~ 1, conf.type = type)
cat(type, ": ", round(min(f$lower), 3), " - ", round(max(f$upper), 3), "\n")
}#> plain : 0 - 1
#> log : 0.124 - 1
#> log-log : 0.06 - 0.676
生存率は0から1の間の値なのに、線形の区間はその外へ出ます。Rは0や1で切って表示しますが、切られた端点は「ここまでの可能性がある」という意味を失っています。
$\log(-\log \hat{S}(t))$変換なら、この問題は起きません。この量は実数全体を動くので、正規近似で区間を作って戻せば必ず0と1の内側に入ります。私なら、小標本の解析や曲線の右端を議論する論文では、conf.typeにlog-logを指定します。イベント数が十分あって推定値が0.2から0.8あたりに収まっているなら、どの方法でもほぼ同じ区間になるので既定のままでかまいません。
ただし、どの方法を使っても漸近正規性に頼っている点は変わりません。リスク集合が数人になった領域では、区間が0と1の内側に収まっていることと、区間が信頼できることは別の話です。
右端で何人が残っているか
実データで見てみます。survivalパッケージのlungは、進行非小細胞肺癌228例の追跡データです。
lu <- survfit(Surv(time, status) ~ 1, data = lung, conf.type = "log-log")
summary(lu, times = c(200, 400, 600, 800, 1000), extend = TRUE)#> time n.risk n.event survival std.err lower 95% CI upper 95% CI
#> 200 144 72 0.6803 0.0311 0.6149 0.737
#> 400 57 54 0.3768 0.0358 0.3069 0.446
#> 600 24 22 0.2136 0.0335 0.1520 0.282
#> 800 8 15 0.0783 0.0246 0.0390 0.135
#> 1000 2 2 0.0503 0.0228 0.0179 0.109
ここは読み違えやすいところです。信頼区間の絶対的な幅だけを見ると、右端のほうが狭くなっています。ところが800日時点の区間は0.039から0.135で、上限は下限の3.5倍です。200日時点は0.615から0.737で1.2倍でした。相対的な不確実性は右へ行くほど大きくなっています。
だから、リスクテーブルなしのKaplan-Meier曲線からは、右端の数値を読み取れません。曲線の見た目は左端も右端も同じ太さの線なので、何人の上に載っているかは線からは分からないからです。ここが、リスクテーブルを添えることが慣行を超えて要件になっている理由です。CONSORT声明も追跡と除外の流れを図で示すことを求めています。
曲線から読めるものと読めないもの
生存期間中央値は、推定値が0.5を下回る最初の時点として読みます。曲線の中ほどで決まるので、リスク集合がまだ大きく、比較的安定した要約値になります。
読めないものが2つあります。1つは最終観察時点より先の生存率です。Kaplan-Meier推定量は観測された時点の上でしか値を持たないので、3年追跡の研究から5年生存率は出せません。曲線が右端で水平に伸びているのは、そこで打ち切りになったという意味であって、その水準が続くという意味ではありません。
もう1つは、競合イベントがあるときの累積発生率です。がん死亡を見たいときに他病死を打ち切り扱いにすると、$1 – \hat{S}(t)$は累積発生率を過大に評価します。この点は打ち切りと左切断で扱いました。
そして、Kaplan-Meier推定量は共変量を調整できません。年齢や病期が群間で偏っている観察研究では、曲線の差をそのまま治療効果と読むことはできません。調整が必要ならCox比例ハザードモデルに進みます。
図と表をどう作るか
投稿用のKaplan-Meier図には、曲線の下にリスク集合数を置きます。目盛りは追跡期間を4から6等分した位置にとり、打ち切りはtick markで示す。ここまでが最低限です。
そのうえで、リスク集合が一定人数を下回る領域をどう扱うかを決めます。曲線を描くのをそこで止める、あるいは描いたうえで本文で言及する。どちらでもかまいませんが、右端の数値を結論に使うなら、そこで何人が残っていたかを本文に書く必要があります。「5年生存率は62%（95% CI 41-78、5年時点のリスク集合7例）」と書いてあれば、読者は自分で重みづけできます。
数値表のほうには、特定時点の生存率を信頼区間つきで載せます。中央値が未到達なら空欄にせず、そう書いたうえで24か月や36か月の生存率で補います。信頼区間の構成方法も脚注に書いておくと、区間の再現性が上がります。
推定量そのものの性質を追うなら、累積ハザードの側から推定するNelson-Aalen推定量が対になります。小標本では両者の値がずれ、Kaplan-Meier推定量のほうが生存確率をやや低めに出す傾向があります。群間比較の検定に進む場合、そして生存曲線が交差している場合に何が起きるかは、比例ハザード仮定の検証で扱います。
自分のデータで、曲線のどこまでを結論の根拠にしてよいか。リスク集合が何人を切ったら本文の書き方を変えるか。ここは追跡の設計とイベント数を見ないと決まりません。Dr.データサイエンスでは、こうした図表の作り方と記載のご相談を承っています。
参考文献
Kaplan EL, Meier P. Nonparametric estimation from incomplete observations. Journal of the American Statistical Association. 1958;53(282):457-481.
Greenwood M. The natural duration of cancer. Reports on Public Health and Medical Subjects. 1926;33:1-26.
Pocock SJ, Clayton TC, Altman DG. Survival plots of time-to-event outcomes in clinical trials: good practice and pitfalls. The Lancet. 2002;359(9318):1686-1689.
Rulli E, Ghilotti F, Biagioli E, et al. Assessment of proportional hazard assumption in aggregate data. British Journal of Cancer. 2018;119(12):1456-1463.
Klein JP, Moeschberger ML. Survival Analysis: Techniques for Censored and Truncated Data. 2nd ed. Springer; 2003.
