RでKaplan-Meier plotを作成する - 統計学入門一歩先へ

はじめに

がんの治療後の生存期間や、心筋梗塞・脳卒中などの疾患の発生までの期間を記述する方法としてよく用いられるのが「Kaplan-Meier plot」と呼ばれるグラフである。例えば下記の図は、乳がんに対する2つの治療法の生存期間をKaplan-Meier plotで要約したものである[1]。

f:id:mstour:20200807195831p:plain — 生存期間を記述するKaplan-Meier plotの例(Im et al.(2019)より抜粋)

このplotの曲線は、「Kaplan-Meier法」によって推定した各時点での「生存関数」 $S(t)$ を描いたものである。生存関数とは、各個人の生存時間を確率変数と考えたときに、その時点まで生存している確率を示す。累積分布関数 $F(t)$ とは逆の関係

$\displaystyle S(t) = 1 - F(t)$

にある。
この手の研究においては、研究の途中で多くの参加者と連絡が取れなくなり、以降の状況が不明になってしまうことが一般的である。このような「打ち切り」があると情報のロスが生じてしまう。Kaplan-Meier法は、打ち切りとなった個人の情報もその時点までの推定に用いて計算を行う。上のplotでは、四角や三角のマークが打ち切りのあった時点を表している。
さて、今回の本題は、Rでこのようなplotを作図する方法である。理論的な部分の説明は、また機会があればさせていただきたい。

Rで作成する。

Kaplan-Meier法による生存関数の推定自体はsurvivalパッケージだけでできるが、survminerパッケージを使うことできれいにplotが描けるのでこちらもあわせて使用することにする（というか、それ以外の方法が思いつかない）。
そういえば、Rを覚えたてのころはパッケージの使い方もいまいちわかっていなくて、library(xxxxx)を書かずにパッケージの関数を使おうとしてエラーが出るのがなぜか分からず延々時間を費やしたり、さらにはインストールしていないパッケージを使おうとしたりしていた・・・
というわけで、パッケージをインストールして（インストール手順は省略）、使う前にちゃんとパッケージを読み込もう（笑）

library( survival )
library( survminer )

さて、まずはKaplan-Meier法による生存関数の推定を行う。これはsurvivalパッケージのsurvfitを使えばよい。

KM = survfit( Surv( TIME, CENSOR ) ~ GROUP, data = DATASET )

ここではKMに推定した生存関数を格納している。結果変数はTIME（生存時間）とCENSOR（打ち切りフラグ：0が打ち切り、1がイベント）をSurv関数で組にしたものを使用し、GROUPで群別したい変数を指定する（ない場合は空白）。DATASETは使用するデータセットである。

次に生存関数をプロットする。ここでsurvminerパッケージのggsurvplotを使用する。

ggsurvplot( KM, data = DATASET )

次のようなplotが作成される。この例では男女別に生存曲線をプロットしている。曲線の途中にいくつも縦棒がくっついているが。これが打ち切りを表している。

f:id:mstour:20200107195538p:plain — Kaplan-Meier plotの例

推定された生存関数は、標本平均などと同様それ自体がばらつきをもった確率変数であるため、信頼区間を構成するこことができる。信頼区間を計算し、プロットに表示するには、以下のようにオプションを追加するだけでよい。

library( survival )
KM = survfit( Surv( TIME, CENSOR ) ~ GROUP, data = DATASET, conf.type = "plain" )

library( survminer )
ggsurvplot( KM, data = DATASET, conf.int = T )

推定量を算出する部分では、survfit関数のconf.typeオプションを使用して、信頼区間を構成する。信頼区間の計算方法はいろいろあって、上の例のように"plain"と指定すると最も古典的な方法であるGreenwoodの公式（デルタ法）で計算した分散にもとづく信頼区間となる（デフォルトがこれだったと思う）。ところがこの方法を使うと信頼区間が0〜1の区間をはみ出す可能性があるとかの問題がある（とは言えよく使われているのを見るのであまり気にしなくていい？）。そのようなことが起こらない方法としてlog-log変換というのもある。この例も下に書いておく。

KM = survfit( Surv( TIME, CENSOR ) ~ GROUP, data = DATASET, conf.type = "log-log" )
ggsurvplot( KM, data = DATASET, conf.int = T )

信頼区間の水準も指定することができるはずだが、95%以外を使うことはあまりないので省略する。

信頼区間のプロットを描くには、ggsurvplotにてconf.int = Tを記載する（デフォルトではF）。以下のようなプロットが作成される。

f:id:mstour:20200108070724p:plain — Kaplan-Meier plotに信頼区間を追加

ついでに細かいことを言うと、このような信頼区間は生存時間（横軸）の各点ごとの生存関数の信頼区間を出していて、生存関数全体が95%の信頼度ですっぽりおさまるような区間ではない（こういうのは信頼帯とか信頼幅とか言う）。

ところで元々はproduct-limit推定量という名前で提案されたようだが、いつからか開発者の名前がそのまま手法の名称になったようだ（イベントの発生した時点における条件付き生存確率の積（product）で生存関数推定量を表し、時間区間を極限（limit）まで小さくすると真の生存関数に収束する、というのが由来だったと思う）。

おわりに

今回はRでKaplan-Meier plotを作成する方法を紹介した。医学研究では知らない人はいないような方法だが、他の分野でどれだけ使われているかは不明である。応用されているような分野があればぜひご教示いただきたい・・・

参考文献

[1] Im S. A. et al.(2019). "Overall Survival with Ribociclib plus Endocrine Therapy in Breast Cancer", N Engl J Med 381:307-316.
[2] デビッド・ホスマー他(2014). "生存時間解析入門原書第2版", 東京大学出版会.