文書の過去の版を表示しています。

目次

ノンパラメトリックな生存関数の推定法

カプラン・マイヤー法は、特定の確率分布を仮定しない生存関数の推定法である。 その推定量(推定値の計算方法)を、カプラン・マイヤー推定量と呼ぶ。 ある時点の生存関数をそれまでの各時点の生存率の積で推定する、という計算方法を用いており、この推定量の性質を示す際にproduct limitという考え方が導入されることから、product limit推定量とも呼ばれる。

survivalパッケージの中の関数survfitは、この推定量を算出することと、グラフを描く情報を提供することができる。survfitを使うには、イベント発生時点とイベントの種類という２種類のフィールド(変数)を持つデータを、データフレームで用意する。

イベント発生時点とは年月日等の暦上の時点ではなく、製品の故障なら使い始めてから故障が起きるまでの時間、人の寿命なら亡くなった時の年齢、治療を開始してから治療を終了するまでの期間など、イベントが起きるまでの期間の長さである。

イベントの種類には故障あるいは打ち切りの２種類がある。故障を1、観測打ち切りを0と、いずれも数値で与えることになっている。製品の故障とは、あらかじめ定められた機能をあらかじめ定められたようには果たせなくなること、人の故障とは、治療中の病気が治る見込みがなく次の段階に進行することや、亡くなること、である。製品の打ち切りとは、その時点で観測を終了してその後の状況を把握できなくなること、人の打ち切りとは、治療のために通っていた病院に何らかの理由や事情で通院しなくなったためにその後の状況を把握できなくなること、である。寿命データ解析では、故障でないことが確認できていた一番最後の時点で、観測を打ち切った、ものとして扱う。

イベント発生時点イベントの種類timestatus
5.5故障5.51
10.3故障10.31
10.0打ち切り5.50
3.5故障5.51

このデータをデータフレームで与えるには、次のように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

推定した生存関数を描くには、関数plotを用いる。 

plot(survfit(Surv(time,status)~1,data=X))

凡例をつけるなら、次の一行を追加で実行すると良い。 

legend("topright",
       lty=c(1,2), 
       legend=c("Kaplan-Meier Estimate", "95% Confidence Limits"))