生存解析とは?
生存解析は興味のある単一のイベントが発生するまでの時間 t に関する分析で、主に2つの分野で利用されています。
一つは医学の分野での利用で、多くの場合治療を開始してからある状況が発生するまで、もしくは患者が死亡するまでの時間に関するものです。 もう一つは工業分野における利用で、製品の信頼性と製品の故障までの時間に関するものです。
nAGでは生存解析に関してカプラン・マイヤー推定量を計算するためのルーチン、コックスの比例ハザードモデルを取り扱う計算ルーチン他を提供しています。 これらの計算ルーチンはExcelをはじめ、MATLAB等のソフトウエアや、C/C++、Fortranなど様々な言語から利用可能です!
生存解析について
故障時間の確率密度関数を f(t) とした場合、最低でも時間 t まで生存する確率を示す生存関数 S(t) は以下の式により与えられます。
上記式で F(t) は累積密度関数を示しています。 ハザード関数 λ(t) はおのおのが t で与えられる時間まで生存した場合に時間 t にて故障が起こる確率を示し、以下の式で与えられます。
累積ハザード率は以下の式で定義されます。
よって
多くの場合、生存解析の対象となるデータは一部が右側打ち切り(例えば故障や死亡などのイベントが発生していない状況)となりますが、このことは解析の計算をより難しくしています。
生存関数の推定
打ち切りデータを含む生存関数の最も一般的な推定はカプラン・マイヤー(Kaplan-Meier)法です。(Product-limit推定)
上記式において dj は、tj で示される時間において nj あるうちのいくつが故障(failure)しているかを示します。これはそれぞれの故障時間(打ち切りは含まれない)のステップを持ったステップ関数です。
であるので、累積ハザード率は以下のように推定できます。
もしくはその対数とtもしくはその対数でプロット図を作成することは、適切なパラメトリックモデルを見つけ出すことに多くの場合役立ちます。 以下に示す関係をその際に利用できます。
- 指数分布:Λ(t)=λt
- ワイブル分布(Weibull Distribution):log λ + γlog t
- ゴンペルツ分布(Gompertz Distribution):log(λ(t))=logλ+γt
- 極値分布:log(Λ(t))=λ(t-γ)
比例ハザードモデル
生存解析を行う多くの場合において、ハザード関数といくつかの説明変数もしくは共変量(covariates)との関係がモデル化されます。 共変量は例えば各々の観察開始時における群もしくは処置の値などがあります。共変量にはこのような時間に依存しないものと時間に依存するものがありますが、時間に依存する共変量は外部共変量と内部共変量に分類されます。この中で外部共変量は故障に直接関係の無いものをそのように呼びます。
ハザード関数に共変量を関連づける最も一般的な関数は比例ハザード関数です。
ここで λ0(t) はベースラインハザード関数で、 z は共変量のベクトルで、 β は未知のパラメータのベクトルです。 また共変量がハザードにたいして複数の影響を及ぼすと言うことが想定されています。
λ0(t) の形式は上記に示された分布のうちの一つもしくはノンパラメトリックな関数であってもかまいません。指数分布、ワイブル分布、極値分布のいずれかである場合には、Aitkin and Clayton (1980)が説明しているポワソン誤差の一般化線形モデルを使った方法で、打ち切りデータへのあてはめ(フィッティング)が行えます。 その他利用可能なモデルにガンマ分布と対数正規分布 があります。
コックスの(Cox's)比例ハザードモデル
特定の形式のハザード関数を使うかわりに、Cox (1972)は λ0(t) が定義されない時間関数の場合を考えました。 このような固定共変量を想定するモデルのあてはめには周辺尤度(marginal likelihood)を使用します。故障が発生したそれぞれの時間 ti のうち ti で打ち切りにされたデータを含む今現在も対象としているデータは ti のリスク集団(risk set)と呼ばれ R(t(i)) で表されます。リスク集団が与えられた場合、故障の可能性があるすべての di の観察対象のうち観察された di 個のケースが故障する確率は以下のように表されます。
ここで si は観察された di 個の個体が t(i) で故障する場合の共変量の和で、 R(t(i)) から抜き出された ni 個の個体から構成されるすべての異なる群を足し合わせたものです。 これは複雑な尤度となりますが、もし故障回数に関連がない場合にはより簡潔に以下のようになります。
ここで nd は異なる故障数の数を表します。 故障回数に関連がある場合には以下の近似が利用可能です。
モデルのあてはめ後に(λ0(t) と残差から)ベースライン生存関数を計算することが可能です。 これを利用してモデルの善し悪し(特に比例ハザードの想定)を調べる事が可能です。
生存解析に関連するnAG数値計算ライブラリのルーチン詳細
G12AAF - カプラン・マイヤ推定値の計算
G12BAF - コックスの比例ハザード・モデルのフィット
G12ZAF - 固定共変量に対するコックスの比例ハザード・モデルに伴うリスク集団の作成
G12ABF - ログランク検定を含む数種の検定量の算出
G01MBF - ミル比の逆数の計算(正規分布におけるハザード率)
G02GCF - ポアソン誤差をもつ一般化線形モデル・フィット
G02GDF - ガンマ誤差をもつ一般化線形モデル・フィット
G07BBF - グループ化データ及び/または打ち切りデータから標準正規分布の母数に対する最大尤度推定値の計算
G07BEF - ワイブル分布の母数に対する最大尤度の推定値の計算
G08RBF - 順位を使った回帰,右打ち切りデータ
G11CAF - 条件付きロジスティックモデル・フィット。(G12ZAFからのリスク集団に利用するとコックス比例ハザードモデルで関連する観察がある場合に周辺尤度を用いたあてはめが行われる)