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

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)
Surv(gehan$time,gehan$cens)

 図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 、共変量を導入しない、分布を仮定しない)
◇セミノンパラメトリックモデル(semi-non-parametric model、共変量を導入する、分布を仮定しない)
◇パラメトリックモデル(parametric model 、共変量を導入する、分布を仮定する)

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

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

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

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

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

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

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

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

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

をネルソン(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
ge.sfの結果

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

>summary(aml.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":
plain


  "log":
log


  "log-lig":
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)
survdiff(Surv(time)~treat,data=gehan)の結果

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