[連載] フリーソフトによるデータ解析・マイニング 第37回

Rと生存時間分析(1)

1.基本概念

(1) 生存時間分析

  生存時間分析 (survival analysis) は、イベント (event) が起きるまでの時間とイベントとの間の関係に焦点を当てる分析方法である。生存時間分析は、工学分野においては機械システムや製品の故障などを、医学分野においては疾患の病気の再発や死亡などを対象とした研究分野である。このような故障、破壊、倒産、再発、死亡などのイベントの生起のことを広義で死亡と呼ぶことにする。

(2) 打ち切り

  生存時間に関する試験・観察を行うとき、治療の中止、もしくは転院、試験・観察の途中で脱落する場合がある。また研究の終了時点では死亡に関するデータが入手できないことが起り得る。このようなことを打ち切りが生じた (censoring) と言う。打ち切りは、左側打ち切り、右側打ち切りなどに分類する場合があるが、ここでは研究を終了するまでイベントが観測できなかったケースを打ち切りと呼ぶことにする。
  Rの生存時間分析のパッケージとして最も広く用いられているのは survival である。パッケージ survival は、Rをインストールする際に自動的にインストールされる。またパッケージ MASS には gehan という生存時間データがある。データ gehan は、白血病患者に対する薬の効果を調べるために被験者42吊に対して行った臨床試験データである。被験者は薬の投与群と対照群(投薬していない群)のペアによって構成されている。生存時間データの形式および打ち切りに関する情報を確認するためデータフレーム gehan の構造を次に示す。

> library(survival);library(MASS)
> data(gehan); dim(gehan);
> gehan[1:6,]
 [1] 42  4

  pair time cens treat
1 1 1 1 control
2 1 10 1 6-MP
3 2 22 1 control
4 2 7 1 6-MP
5 3 3 1 control
6 3 32 0 6-MP

  変数 pair は投薬と比較対象のペア、time は生存時間、cens は打ち切りか否か、 treat は 6-MP 薬の投与か否かである。打ち切りか否かは0,1で表している。
  パッケージ survival には、生存時間データオブジェクトを作成する関数 Surv があり、返された結果は通常生存時間モデルの目的変数(あるいは被説明変数)として用いられる。
  次に、関数 Surv を用いた例を示す。返される結果の値の右側に記号+が付いている値は打ち切りデータである。

> Surv(gehan$time,gehan$cens)

[1] 1 10 22  7  3 32+ 12 23  8 22 17  6  2 16 11 34+ 8
[18] 32+ 12 25+ 2 11+ 5 20+ 4 19+ 15  6  8 17+ 23 35+ 5  6
[35] 11 13  4  9+ 1  6+ 8 10+

  図1に横軸を生存時間としたデータ gehan の棒グラフを示す。棒の右側の×印は、研究が終了する前に死亡が確認されたもので、△印は研究が終了する前に何らかの理由で観測が継続できなかった打ち切りであり、○印は研究が終了する時点まで生存した打ち切りである。打ち切りがあるデータを上完全データ、打ち切りがないデータを完全データと呼ぶ。

図1 生存時間の横棒グラフ  

図1 生存時間の横棒グラフ

(3) 生存関数とハザード関数

  生存時間T を確率変数、確率密度関数を

確率密度関数
累積確率分布関数をF (t ) で表すと、イベントがある時点t まで生起していない生存関数S (t )(survival function)は

生存関数S(t)

で表わされる。
  イベントがある時点t までに生起していないという条件の下で,次の瞬間にイベントが生起する瞬間死亡率を示すハザード関数h (t ) (hazard function) は

ハザード関数

となる。ハザードは危険度とも呼ばれている。生存関数S (t ) とハザード関数h (t ) および累積ハザード関数H (t ) との関係は、

生存関数とハザード関数および累積ハザード関数との関係

のように表わすことができる。

(4) 生存時間分析の分類

  生存時間分析は、生存時間に影響を与える時間以外の共変量(複数の要因、説明変数)がパラメータとして作成するモデルに導入されているか否か,生存時間の分布形に特定の確率分布を仮定するか否かによって、次のように分類することができる。

  @ ノンパラメトリックモデル (non-parametric model 、共変量を導入しない、分布を仮定しない)
  A セミノンパラメトリックモデル (semi-non-parametric model、共変量を導入する、分布を仮定しない)
  B パラメトリックモデル (parametric model 、共変量を導入する、分布を仮定する)

2.ノンパラメトリックモデル

(1) ノンパラメトリックの推定法

  ノンパラメトリックの推定方法とは、確率分布を仮定せずに生存時間を推定する方法で、経験分布による推定法とハザード関数による推定法がある。

  @ 経験分布による推定法
  打ち切りがない完全データであれば、収集したデータに基づいた累積分布関数、つまり経験分布関数F (t ) を用いて推定することができる。

平均S(t)=1-F(t)

  これは後述のカプラン・マイヤー (Kaplan Meier) 推定法の特殊なケースである。
  カプラン・マイヤー推定法は、条件付き確率の考え方に基づき、次のように定義される。

カプラン・マイヤー推定法の定義

  式の中の d i は、時点 t i の死亡数、r i は時点 t i の直前までイベントが生じる可能性のある観察対象者の数 (リスクセット) である。

  A ハザード関数による推定法
  カプラン・マイヤー推定量を用いて累積ハザード関数を推定することができる。累積ハザード関数およびその推定量はそれぞれ

累積ハザード関数およびその推定量

である。カプラン・マイヤー推定量から求めたハザードの推定量

ハザードの推定量

をネルソン (Nelson) 推定量、あるいはネルソン・アーラン (Nelson-Aalen) 推定量と呼ぶ。これは小サンプルに有効であると言われている。ネルソン推定量を修正した

フレミング・ハリントン(Fleming-Harrington)推定量

をフレミング・ハリントン (Fleming-Harrington) 推定量と呼ぶ。

(2) 関数 survfit

  パッケージ survival には、ノンパラメトリック法による生存時間を当てはめる関数 survfit がある。関数 survfit は、引数 type を用いて推定法を指定する。関数 survfit には3種類の推定方法が用意されている。デフォルトに設定されているカプラン・マイヤー推定法 (Kaplan*Meier) のほかに Fleming*Harrington 推定法、fh2 推定法がある。関数 survfit の書き式を次に示す。

survfit(formula, data,type=" ",…)

  引数 formula では、次に示す例のように Surv オブジェクト形式の目的変数と説明変数を指定する。次にデータ gehan の生存時間 time を目的変数、投薬したか否かのデータを記録した treat 列を説明変数とした関数 survfit の使用例を示す。その他の引数はデフォルトの値を用いた。

> ge.sf<-survfit(Surv(time,cens)~treat, data=gehan)
> ge.sf
Call: survfit(formula = Surv(time, cens) ~ treat, data = gehan)

n events median 0.95LCL 0.95UCL
treat=6-MP
21 9 23 16 NA
treat=control
21 21 8 4 12

  このように、データ gehan における投薬群と対照群ごとの標本数、イベントの件数、中央値、両側の区間推定に関する情報を返す。
  関数 summary を用いると、各標本の生存時間 time、リスクセット n.risk、イベントの数 n.event、推定された生存確率 survival、標準誤差 std.err、95%信頼区間の上下限値を返す。

> summary(ge.sf)
<結果は省略>

  次のように関数 plot を用いると生存確率の推測値のプロットを返す。

> plot(ge.sf,lty=1:2)
> legend(locator(1),c("6-MP投与群","対照群"), lty=c(1,2))

図2 データgehanの生存時間プロット

図2 データ gehan の生存時間プロット

  生存曲線プロットの折れ線に+印が付いているところは打ち切りがある時点である。関数 plot に引数 mark.t=F を用いると印+が除かれる。また、引数 formula で、打ち切りデータを省略することができる。作図する時、引数 conf.int = TRUE を加えると生存曲線プロットに信頼区間が表示される。ただし、生存曲線が1本である場合は、conf.int = TRUE を省略しても、信頼区間が作成される。
  次に gehan の投薬群のみのデータについてカプラン・マイヤー法による生存曲線と90%の信頼区間を示す作図コマンドの例を示し、その結果を図3に示す。

> ge2<-subset(gehan,treat=="6-MP")
> ge2.s <-survfit(Surv(time,cens)~treat, conf.int=.9,data= ge2)
> plot(ge2.s,mark.t=F)
> legend(locator(1),lty=c(1,2),legend=c("生存曲線","90%信頼区間"))

図3 生存曲線と信頼区間

図3 生存曲線と信頼区間

  カプラン・マイヤー推定量 は漸近的に正規分布に従うことが知られており、標準偏差は、次の式を用いて推定することができる。

標準偏差の推定

  ネルソン推定量の標準偏差は次の式で求める。

ネルソン推定量の標準偏差

  関数 survfit には、信頼区間を推定する方法を指定する引数 conf.type が用意されている。引数 conf.type は、"none", "plain", "log", "log-log" の中から1つ選択できる。デフォルトには"log"になっている。"none"の場合は、信頼区間を返さない。それぞれの conf.type は次の式で推測される。式の中の S は生存関数の推定値 S(t)H はハザード関数の推定値H(t)se は標準誤差を示す。

  plain,log,log-log

  "plain" と "log" で求めた信頼区間の上限値は1を超える場合があるが、そのときには1を超える部分は切り捨てる。"log-log" による計算結果は1以内に収まる。
  データ ge2 の3種類 ("plain", "log", "log-log") の conf.type の信頼区間を求め、図示するコマンドの例を次に示し、その結果を図4に示す。

> plot(ge2.s,conf.int=TRUE,mark.t=F)
> lines(survfit(Surv(time,cens)~treat, data=ge2,conf.type="plain"),mark.t=F, conf.int=TRUE,lty=3,col=3)
> lines(survfit(Surv(time,cens)~treat, data=ge2,conf.type="log-log"),mark.t=F, conf.int=TRUE,lty=4,col=4)
> legend(locator(1),c("log","plain","log-log"), lty=c(1,3,4),col=c(1,3,4))

図4 異なる推定法による信頼区間

図4 異なる推定法による信頼区間

  関数 survfit には、信頼区間を指定する引数 conf.int があり、デフォルトは conf.int = .95 になっている (95%の信頼区間)。引数 conf.int の値は自由に指定することができる。
  関数 survfit における引数 type に fleming-harrington (fh に略することができる) あるいは fh2 を指定するとフレミング・ハリントン推定量を返す。異なる推定法の結果を考察する際に必要となる生存曲線のプロットを作成するコマンドの例を次に示し、その結果を図5に示す。

> ge2.s <-survfit(Surv(time,cens)~ treat, data= ge2)
> ge2.fh<-survfit(Surv(time,cens)~treat, data= ge2,type="fleming-harrington")
> ge2.fh2<- survfit(Surv(time,cens)~treat, data= ge2,type="fh2")
> plot(ge2.s,conf.int=F,mark.t=F)
> lines(ge2.fh,lty=2)
> lines(ge2.fh2,lty=3,col=2)
> legend(locator(1),lty=1:3,legend=c("Kaplan-meier","fleming-harrington","fh2")

図5 3種類の推定法による生存曲線

図5 3種類の推定法による生存曲線

(3) 生存関数の検定

  データ gehan のような2群 (投薬群、対照群など) 以上の観測値が得られたときには、群ごとの生存曲線の間の差の有意性について検定を必要とする場合がある。最も広く用いられている検定方法はログ・ランク (log-rank) 検定法である。ログ・ランク (log-rank) 検定法は、群ごとのイベントありとなしに集計したクロス表のカイ2乗検値の推定値を検定統計量とする。
  パッケージ survival には、検定を行う関数 survdiff がある。関数 survdiff は、G-rho ファミリの検定を行う。引数 rho = 0 にするとログ・ランク検定、rho = 1 にすると Gehan-Wilcoxon 検定を行う。デフォルトは rho = 0 になっている。次にデータ gehan の生存曲線のログ・ランク検定の例を示す。

> survdiff(Surv(time)~treat,data=gehan)

Call:
survdiff(formula = Surv(time) ~ treat, data = gehan)

N Observed Expected (O-E)^2/E (O-E)^2/V
treat=6-MP
21 21 29.2 2.31 8.97
treat=control
21 21 12.8 5.27 8.97

Chisq= 9 on 1 degrees of freedom, p= 0.00275

  ログ・ランク検定の p 値は約0.003であるので、仮に通常よく用いられている有意水準5%を規準とすると、両群 (投薬群と対照群) の生存曲線には有意な差が認められる。
  次の号では、セミノンパラメトリックモデル、パラメトリックモデルについて紹介する。