カプラン・マイヤー法
生存関数は時点ごとの生存率を表す関数であり、カプラン・マイヤー法は、特定の確率分布を仮定しない生存関数の推定法である。ノンパラメトリックな生存関数の推定法と言う。 その推定量(推定値の計算方法)を、カプラン・マイヤー推定量と呼ぶ。 ある時点の生存関数をそれまでの各時点の生存率の積で推定する、という計算方法を用いており、この推定量の性質を示す際にproduct limitという考え方が導入されることから、product limit推定量とも呼ばれる。
survivalパッケージの中の関数survfitは、この推定量を算出することと、グラフを描く情報を提供することができる。survfitを使うには、イベント発生時点とイベントの種類という２種類のフィールド(変数)を持つデータを、データフレームで用意する。
イベント発生時点とは年月日等の暦上の時点ではなく、製品の故障なら使い始めてから故障が起きるまでの時間、人の寿命なら亡くなった時の年齢、治療を開始してから治療を終了するまでの期間など、イベントが起きるまでの期間の長さである。
イベントの種類には故障あるいは打ち切りの２種類がある。故障を1、観測打ち切りを0と、いずれも数値で与えることになっている。製品の故障とは、あらかじめ定められた機能をあらかじめ定められたようには果たせなくなること、人の故障とは、治療中の病気が治る見込みがなく次の段階に進行することや、亡くなること、である。製品の打ち切りとは、その時点で観測を終了してその後の状況を把握できなくなること、人の打ち切りとは、治療のために通っていた病院に何らかの理由や事情で通院しなくなったためにその後の状況を把握できなくなること、である。寿命データ解析では、故障でないことが確認できていた一番最後の時点で、観測を打ち切った、ものとして扱う。
| イベント発生時点 | イベントの種類 | time | status | |
| 5.5 | 故障 | 5.5 | 1 | |
| 10.3 | 故障 | 10.3 | 1 | |
| 10.0 | 打ち切り | 10.0 | 0 | |
| 3.5 | 故障 | 5.5 | 1 | |
このデータを小さい順に並び替える。
| イベント発生時点 | イベントの種類 |
| 3.5 | 故障 |
| 5.5 | 故障 |
| 10.0 | 打ち切り |
| 10.3 | 故障 |
そして時点ごとに、それまでに生存していた数、故障数、打ち切り数を求める。
| イベント発生時点 | イベントの種類 | 時点直前までの生存数 | その時点での故障数 | その時点での打ち切り数 |
| 3.5 | 故障 | 4 | 1 | 0 |
| 5.5 | 故障 | 3 | 1 | 0 |
| 10.0 | 打ち切り | 2 | 0 | 1 |
| 10.3 | 故障 | 1 | 1 | 0 |
各時点での故障率は、故障数/生存数で推定できる。
| イベント発生時点 | イベントの種類 | 時点直前までの生存数 | その時点での故障数 | その時点での打ち切り数 | 故障率 |
| 3.5 | 故障 | 4 | 1 | 0 | 1/4 |
| 5.5 | 故障 | 3 | 1 | 0 | 1/3 |
| 10.0 | 打ち切り | 2 | 0 | 1 | 0/2 |
| 10.3 | 故障 | 1 | 1 | 0 | 1/1 |
生存率は、前の時点の生存率×(1-その時点の故障率)で繰り返し計算できる。
| イベント発生時点 | イベントの種類 | 時点直前までの生存数 | その時点での故障数 | その時点での打ち切り数 | 故障率 | 生存率 |
| 3.5 | 故障 | 4 | 1 | 0 | 1/4 | 1×(1-1/4)=0.75 |
| 5.5 | 故障 | 3 | 1 | 0 | 1/3 | 0.75×(1-1/3)=0.50 |
| 10.0 | 打ち切り | 2 | 0 | 1 | 0/2 | 0.50×(1-0/2)=0.50 |
| 10.3 | 故障 | 1 | 1 | 0 | 1/1 | 0.50×(1-1/1)=0 |
これがカプラン・マイヤー法である。 簡単な計算のように見えて、この方法の性質を検討するには、少し難しい確率論を学ぶことになる。 以下では、この方法が正しいとして、話を進める。
関数survfit
このデータをデータフレームで与えるには、次のようにRに伝える。
X = data.frame(time=c(5.5, 10.3, 10.0, 3.5), status=c(1, 1, 0, 1))
日本語の変数名を使える設定にしてあれば、次のようにしてもいい。
X = data.frame(時間=c(5.5, 10.3, 10.0, 3.5), イベント=c(1, 1, 0, 1))
このデータからカプラン・マイヤー法で生存関数を推定するには、関数survfitを次のように用いる。
survfit(Surv(time, status)~1, data=X)
survfitの一つ目の引数は、モデル式と呼ばれる。データフレームはdata=という引数で与える。その他、survfitの詳細は次のとおり。
| 関数 | survfit |
| パッケージ | survival |
| 引数 | モデル式, データフレーム名 |
| モデル式の書き方 | Surv(時間, イベント種類)~1 |
| データフレーム名の与え方 | data= |
| 返り値(表示) | Call:survfitが実行されたときのコマンド, n:サンプル数, events:故障数, median:寿命のメディアンの推定値, 0.95LCL:メディアンの下側95%信頼限界, 0.95UCL:メディアンの上側95%信頼限界 |
| 返り値(詳細) | n:サンプル数, time:イベント発生時点, n.risk:各時点のリスクセットの大きさ(故障率の分母), n.event:各時点の故障数(故障率の分子), n.censor:各時点の打ち切り数(故障率の分母の減少分), surv:各時点の生存関数の推定値, type:打ち切りの種類, std.err:各時点の生存関数の標準誤差, lower:各時点の生存関数の下側信頼限界, upper:各時点の生存関数の上側信頼限界, conf.type:信頼区間の作り方, conf.int:信頼係数, call:survfitが実行されたときのコマンド |
実行すると次の表示を得る。
Call: survfit(formula = Surv(time, status) ~ 1, data = X)
n events median 0.95LCL 0.95UCL
4.0 3.0 7.9 3.5 NA
より詳細な情報は、関数summaryを用いて抽出できる。
survfit(formula = Surv(time, status) ~ 1, data = X)
次のように、データが各イベントの時点とそこでの生存関数を計算するための表に整理して表示される。
Call: survfit(formula = Surv(time, status) ~ 1, data = X) time n.risk n.event survival std.err lower 95% CI upper 95% CI 3.5 4 1 0.75 0.217 0.426 1 5.5 3 1 0.50 0.250 0.188 1 10.3 1 1 0.00 NaN NA NA
この表示を読みながら、カプラン・マイヤー法の説明を少し加える。
- データではイベント発生時点は小さい順に3.5, 5.5, 10.0, 10.3の4時点だった。それが上の表示では10.0がなくなっている。
- 時点3.5までは4つの対象が生存している、あるいは故障していない。時点3.5で一つ故障したので、それ以降、次の時点5.5までは3つの対象が生存している、あるいは故障していない。これがn.riskという列の数値の意味である。時点10.3で1とあるのは、それ以前の時点10.0で対象の1つの観測が打ち切られたためである。時点10.3の故障発生時に、それまで生存していることが確認できているのは1つだけだから。
- 各時点の故障率は1/4、1/3、1/1である。これはn.event/n.riskという割り算で求まる。
- survivalは生存関数の推定値である。時点3.5より前には、すべてが生存していたので、時点3.5の故障によって生存率が1から1×(1-1/4)=0.75に減少した。次の時点5.5では、更に0.75から0.75×*(1-1/3)=0.25に生存率が減少した。このように計算するのが、カプラン・マイヤー法である。
- std.errはグリーンウッドの公式など、指定された方法を用いている。信頼区間も同様で、ここでは詳細を省略する。
推定した生存関数を描くには、関数plotを用いる。
plot(survfit(Surv(time,status)~1,data=X))
凡例をつけるなら、次の一行を追加で実行すると良い。
legend("topright",
lty=c(1,2),
legend=c("Kaplan-Meier Estimate", "95% Confidence Limits"))
実線の折れ線が、カプラン・マイヤー法で推定した生存関数である。1から生存関数を引けば、累積分布関数の推定ができる。このグラフからまずは、寿命分布の密度関数が一山か、二山以上か、を読み取ろうとする。
survfitが出力した各時点でのリスクセットの大きさn.riskとその時点での故障数n.eventの比を計算すると、故障率曲線も推定できる。
X.survfit = survfit(Surv(time, status)~1, data=X) X.survfit$n.event/X.survfit$n.risk
これを縦軸に、時間を横軸にとったグラフ
plot(X.survfit$time,X.survfit$n.event/X.survfit$n.risk,
xlab="Time", ylab="Failure Rate")
が故障率曲線のノンパラメトリックな推定値となる。
ノンパラメトリックな推定量の長所と短所
ノンパラメトリックな推定量は、確率分布を仮定しないので、当初から色眼鏡で見ることがなく、データをそのまま吟味できるという長所を持つ。実際に観測された時点から先の予測ができない、という短所も合わせ持つ。また、精度よく確率分布を推定するには、とても大きいサンプルを必要とすることも、短所である。