大阪電気通信大学

研究論集（自然科学編）第 48 号

フリーソフト R を活用した生存データ解析 ( 総合報告 )Survival Data Analysis Using R
辻谷将明 *1, 田中祐輔 *2
The Cox’s proportional hazards model has been widely used for the analysis of treatment andprognostic effects with censored survival data. The model was developed based on the hazard
function for an individual as a fixed-time covariates. When the covariates values change for theduration of the study, however, a number of theoretical problems that are to be solved with respect
to the baseline survival function and the baseline cumulative hazard function, are involved. This
article presents time-dependent survival data analysis using R which may be evaluated for theirimpact on survival. We consider how these analyses can affect the prediction of patient outcomeusing multistate model, stratified Cox’ proportional hazard model, and competing risk model.
Key words:Cox’s proportional hazard model, Event-history analysis, competing risk, stratified
proportional hazard model, cumulative incidence function
キーワード :Cox 比例ハザードモデル , イベント‐ヒストリー解析 , 競合リスク , 層別比例ハザードモデ
ル , 累積発生関数

１ はじめに
従来、生存時間解析には Cox 比例ハザードモデルが広範に活用されてきた ( 大橋 , 浜田 ,1995; 中村 ,2001;Kalbﬂeish and Prentice, 2002;Klein and Moeschberger, 2003;Lee et al.,2003)。Cox 比例ハザードモデルでは、生存時間 ( 観測時点 ) のハザード関数をと表す。ここで、ベースラインハザード関数

は観測時点の関数であって、共変量を含んでいない。このモデルでは、共変量の観測は各患者について 1 回のみで、時間が変化しても不変である。しかし、反復して測定される検査値の生存時間への影響の評価、あるいは観測期間中に薬剤の投与量を変化させたとき、その効果の有無の検証が必要な場面もでてくる (Andersen et al., 1982;
Aitkin et al., 1983; Altman et al., 1984; Christensen et al., 1986; Aalen et al., 2008)。
本稿では近年 、広範に活用されているフリー・ソフト R( 辻谷 , 外山 , 2007; 辻谷 , 竹澤 , 2009; 辻谷 ,和田 , 2012) を援用し、共変量が時間と共に変動する“時間依存型”の典型的な生存データについて考察する。そして、それを拡張した multi-state モデルに基づくイベントヒストリー解析を取上げる。同一の個体について、イベントが複数回、繰返して発生する場合である。例えば、膀胱癌では再発を繰返すことがよくある。骨髄移植のように移植→再発→死亡と多段階の state を繰返すこともある。次に、競合リスクモデルについて解説する ( 西川 ,2008)。ある個体について、注目している死亡が発生するまでの観測過程で、それ以外の死亡 ( これを競合リスクと呼ぶ ) が起こることがある。これは、従来の比例ハザードモデルを採用し、データを時間依存型にして解析できる。累積発生関数に基づく２群間の比*1*2

大阪電気通信大学 情報通信工学部 情報工学科 , 連絡先〒 572-8530 寝屋川市初町 18-8;E-mail:ekaaf900@ricv.zaq.ne.jp
イーピエス株式会社臨床情報本部データサイエンスセンター , 〒 112-0004 大阪市淀川区宮原 3-4-30 ニッセイ新大阪ビル11 階
2013 年 4 月 30 日受理― 49 ―

This document is provided by JAXA.

較や有意性検定を行う。更に、共変量を考慮し、部分分布とういう考えを取込んだ競合リスク回帰モデ
ルを取上げる。ここでも、共変量に時間依存型データの含まれる場面に拡張できる。

２ 比例ハザードモデル
生存時間を確率変数としたとき，時間 まで生存する確率が，生存時間関数である．累積率分布関数はとなる． の直前まで生存した人が，次のΔ の期間に死亡する条件付き確率は，Δ に依存する．このとき，単位時間当たりの平均死亡率について，のとき，時間 におけるハザード関数をと定義する．すなわち， まで生存した人のうち，に換算し，までに死ぬ人の割合を，単位時間当たりの量としたときの極限値である．ハザード関数

確率密度関数との間にはという関係がある．
個体が生存時間に影響を与える因子 ( 予測変数，背景因子，予後因子 ) を
この値は時間に依存しない．このとき，予測変数と書く．ここに，の予測変数

より，をもつ個体のハザード関数

をベースラインハザード関数，とする．
をを相対リスク (relative risk) と呼ぶ．2 つについてとは比例する ( 比例ハザード性 )．更に，対数線形性を仮定すると比例ハザードモデル― 50 ―

This document is provided by JAXA.

となる．特にならとなる．
(11) 式の未知パラメータは部分尤度を最大化して求める (Cox, 1975)．ここに，

で、は時点でのリスク集合 ( 時点 の直前まで生存した個体から成る集合 ) である．
表 1 は原発性胆汁性肝硬変 (PBC) データである (Collett,2003)。共変量としてビリルビンの初診時の対数をとる。
表 1 原発性胆汁性肝硬変 (PBC) データ

この部分尤度の計算は、表 2 のようになる。
表 2 部分尤度の計算表

通常の Cox の比例ハザードモデルのプログラムは、― 51 ―

This document is provided by JAXA.

となり、解析結果を得る。ビリルビン ( の対数値 ) は、10 ％で有意である。

３ 時間依存型モデルとイベントヒストリー解析
３．
１ 時間依存型モデル
表 1 において、例えば、患者＃ 7 は 1514 日目で死亡するまで、７回来院 (clinic visit) しており、ビリルビン ( の対数 ) 値は、表 3 のように時間と共に変動している。本稿では、このように個体が複数個の観測値を生成している場合を“時間依存型”と呼ぶ。近年、共変量の非線形性を考慮したニューラルネットワークモデル (Tsujitani and Koshimizu,2000;Tsujitani and Sakon,2009,Tsujitani,Iba,Tanaka,2012)、サポートベクターマシン (Tsujitani and Tanaka,2011)、一般化加法モデル (Tsujitaniand Baesens,2012;Tsujitani,Tanaka and Sakon,2012; Tsujitani and Tanaka, 2011; wood, 2000,2006, 2008) などによる解析も考案されている。
表 3 患者＃ 7 に関する時間依存型共変量全ての患者について、最終生存時間の短い順に並べると表 4 のようになる。来院ごとの ln(bil) を示している。最初の死亡例 #281、2 番目の死亡例 #384 に対する部分尤度の計算は表 5 のようになる。
最初の死亡例の部分尤度は― 52 ―

This document is provided by JAXA.

となる。
表 4 全患者に関する時間依存型共変量

表 5 最初の死亡例＃ 281 および 2 番目の死亡例＃ 383 に対する部分尤度の計算

一般的に尤度はから計算される (Lawless,2003,7.1.8 節 )。ここに、は、すべての共変量の値である。

R プログラムは― 53 ―

This document is provided by JAXA.

となり、解析結果を得る。ビリルビン ( の対数値 ) は、尤度比検定 (Likelihood ration test)、およびスコア検定 (Scoretest) で 1 ％有意となる。

３．
２ イベントヒストリー解析
時間依存型の特殊なケースとして、multistate モデルに基づくイベントヒストリー解析がある ( 赤澤ら，2010，第 9 章 )。最も単純な場面は、図 1 のように表せる。観測開始から目標事象 ( 死亡あるいは打切り ) 発生までに種々のイベント ( すなわち，multistate) が起こり、その時間と共変量の値が記録される。

図 1 multistate モデル

膀胱癌データについて、患者 #15，43，86 に関する実際のデータを表 6 に示す。
表 6 患者 #15，43，86 に関する実際のデータ(1)Anderson-Gill モデル
Anderson-Gill モデルでは、表 7 のような時間依存型で表す。Num,Size,Treat は時間と共に変動し― 54 ―

This document is provided by JAXA.

ても良い。例えば、患者 #43 は、時間 3,15,46,51 で再発し、観測終了の 53 時間で打切りとなった。
表 7 Anderson-Gill モデル

R プログラムはとなり、解析結果を得る。よって、Num と Treat が高度に有意となる。
(2)PWP モデル
PWP (Prentice,Williams, and Peterson,1981) モデルは、層別解析の援用である。来院ごとに層別し、データを表 8 のように並び替える。

― 55 ―

This document is provided by JAXA.

表 8 PWP モデル

同表から分るように、来院ごとに層別されており、その層の中では、観測値は独立になり、対数尤度をと分解できる。
R プログラムはとなり、解析結果― 56 ―

This document is provided by JAXA.

を得る。
(3)PWP gap モデル
ＰＷＰモデルで生存時間を時間依存型で与えたが、それをギャップ time=stop-start で置換えたモデルである。R プログラムはとなり、解析結果を得る。

４ 競合リスクモデル
例えば、喫煙習慣と肺癌の疫学調査において、死因として肺癌，その他の癌，その他の疾患，事故が考えられる。このとき、肺癌がエンドポイントで、その他の癌，その他の疾患，事故を競合リスクという。
4 つの死因は ( 確率的に ) 独立ではない。すなわち、肺癌のリスクの高い人は、他の癌のリスクも高い。
打切りのタイプとして
1) 研究 ( 治験 ) の終了：打切りと死亡は独立
2) 追跡不能3) 競合リスク：研究対象以外の死因が発生がある (Putter et al., 2007)。
表 9 は、マウス発癌性データ ( 単位：日数 ) である (Kalbﬂeisch et al., 2002,pp.257-259)。死因とし― 57 ―

This document is provided by JAXA.

て、胸腺リンパ腫、網状組織細胞肉腫、その他がある。共変量は 2 値 ( 対照群と，無菌群 ) である。
表 9 マウス発癌性データ

４．
１ 周辺モデル
表 9 において、周辺モデルでは胸腺リンパ腫を注目しているイベントとし、他の２つの死因を打切りとして比例ハザードモデルを適用する。R プログラムはとなり、解析結果を得る。3 つの死因別の結果を表 10 に与えておく。注目しているイベントを胸腺リンパ腫とした場合のみ、対照群と無菌群の有意差はない。
表 10 死因別の結果

これら 3 種類のイベントごとの累積ハザードを求めると図 2 ～図 4 のようになる。

― 58 ―

This document is provided by JAXA.

図 2 胸腺リンパ腫

図 3 網状組織細胞肉腫

図 4 その他― 59 ―

This document is provided by JAXA.

４．
２ 層別比例ハザードモデル
表 11 は、骨髄腫症データ (Allison,1995,p.31) である。
で、共変量として薬剤に関するは生存時間、を取上げる。

表 11 骨髄腫症データ

生存時間 (dur) と打切り (status) について、比例ハザードモデルを当てはめるととなり、共変量 ( 薬剤 ) は有意にならない。
しかし、患者に陣機能障害があるか否かによってハザード関数の形状が異なることもある。そこで、
層別比例ハザードモデルを採用する。

個の層がある場合、一般にと書け、ベースラインハザード関数は( 層 ) ごとに異なるが、同一の相対リスク関数をもつ比例ハザードモデルである。
層ごとに求めた部分尤度の積の対数― 60 ―

This document is provided by JAXA.

を最大にするが推定量である。データは、表 12 のように与えられる。ここに、とする。部分尤度の計算は表 13 で与えられる。
表 12 腎機能障害で層別された骨髄腫症データ

表 13 部分尤度の計算

R プログラムを用いると― 61 ―

This document is provided by JAXA.

が得られる。層別比例ハザードモデルを採用すると薬剤は有意になる。
次に、共変量が連続値である表 14 のデータを解析する。腎障害，骨病変，タンパク尿などを伴わない良性単クローン性グロブリン血症 (MGUS) に関するデータで、一部が多発性骨髄腫に移行するため，定期的な経過観察が必要である (Therneau et. al., 2000,8.4.1 節 )。共変量として、性別、年齢、ヘモグロビンレベル (hgb)、単クローン性蛋白ピーク (mspike) をもつ。死因としては、MGUS, 骨髄腫、その他の 3 つがある。この死因 ( 競合リスク ) を層別因子とし、表 15 のように書き換える。これは、と書ける。通常の層別比例ハザードモデルとの相違点は、患者が必ず層別されたすべての対数尤度に入ることである ( 通常の層別比例ハザードモデルでは、患者を層ごとに分類する )。
表 14 良性単クローン性グロブリン血症

表 15 競合リスクを層別因子

R プログラムは、となり、解析結果― 62 ―

This document is provided by JAXA.

を得、死因 ( 競合リスク ) をすべてイベント (status=1) として解析した結果と同一になる。

４．
３ 累積発生関数(1) ハザード関数
個の競合リスクの中で、タイプ(Cumulativeのイベントに対する累積発生関数

Incidence Function) はと定義される。すなわち、時間 t 以前に、タイプのイベントが発生する同時確率である。累積発生関数は、部分分布 (subdistribution) とも呼ばれている。競合リスクイベントは互いに独立である必要はない。
時間 以前に、いずれかのイベントが起こる確率は、(20) 式からとなる。時間 までにタイプ のイベントが起こらない確率 ( 部分生存関数という ) は、(20) 式からで与えられる。タイプ のイベントに対する部分密度 (subdensity) 関数は、(20) 式からとなる。
(23) 式から、死因別 (cause-specific) ハザード関数を― 63 ―

This document is provided by JAXA.

と定義する。(24) 式からとの間には (7) 式と類似の関係とがあることに留意されたい。死因別ハザード関数に基づくモデル

が、4.1 節の周辺モデルである。(24) 式からと定義する。(20),(25) 式からとなり、とを推定すれば、が得られる。ちなみにであることに留意されたい。
全死因での確率密度関数からを得る。

より、は分布関数にならない。そのため部分分布と呼ぶ。

部分分布のハザード関数をと定義する (Gray,1988;Pintilie,2006,p.46)。注目しているイベントが、その時間までに起こらないか、あるいは競合リスクイベントが観測された条件のもとで、注目しているイベントが、次の時間こる確率を単位時間当りの量に換算し、に起としたときの極限値が部分分布のハザード関数である。

競合リスクが無ければ ( 注目しているイベントのみ )、死因別ハザードと部分分布のハザードは同一に― 64 ―

This document is provided by JAXA.

なる。
(2) 死因別ハザード関数に基づく CIF のノンパラメトリック推定
イベントが発生した時間をイベント数をとしたとき、時間でのリスク集合の個体数をKaplan-Meier 推定量をとする。

で起きたタイプの、時間 までにいかなるイベントも起きない確率のの直前までにはいかなるイベントも起こらず、 でタイプのイベントを経験する同時確率である CIF はから推定される (Pintilie,2006,4.2.1 節 )。ここに、は、時間でのイベント に対する死因別ハザード関数の推定量である。
区間で個の打切りがで起きたとき、尤度関数はで与えられる (Kalbﬂeisch&Prentice,2002,(8.9),(8.10) 式 )。ただし、プのイベントに対するハザード関数で、は時間における任意のタイとする。(34) 式からを得る。よってよりを得、が求まる。となる。

より、(36) 式はより、(36)、(37) 式を用いるとを得る。(33) 式へ (38) 式を代入すると、は― 65 ―

This document is provided by JAXA.

から推定される。ここに、は、すべてのイベントを同じタイプと見なして得られる Kaplan-Meier推定量である。
(39) 式からを得る (Aalen,1978;Pintilie,2006,4.2.4 節 )。よって、

で与えられる。ここに、は標準正規分布の信頼区間は点である。

計算手順を示すため、数値例として表 16 のデータを取上げる。ただし、イベントのタイプは、0= 打切り、1= 注目しているイベント、2= 競合リスクイベントとする。に算出され、は、表 17 のようとなることを確認できる。
表 16 数値例

表 17 計算手順― 66 ―

This document is provided by JAXA.

R プログラム (Gray,2010) はと書け、図 5 の CIF が得られる。

図 5 CIF の推定

ちなみに、競合リスクを従来のように打切りとして Kaplan-Meier 推定を行うと表 18 のようになる。
図 5 に表 18 のを加えると図 6 のようになる。同図から競合リスクを従来のように打切りとした K-M 推定はナイーブになっていることが分かる。
表 18 競合リスクを打切りとした Kaplan-Meier 推定― 67 ―

This document is provided by JAXA.

図 6 図 5 に表 18 のを加えた CIF

次に、表 9 のマウス発癌性データを解析する。R プログラムはと書ける。イベントごとの CIF は図 7 ～ 図 9 のようになる。

― 68 ―

This document is provided by JAXA.

図 7 胸腺リンパ腫

図 8 網状組織細胞肉腫

図 9 その他(3)CIF に関する 2 群間の差の Gray 検定
Gray(1988) は CIF に関する 2 群間の差を検定するため、部分分布ハザードに関する重み付き平均― 69 ―

This document is provided by JAXA.

を採用した。ここに、とする。
2 群での生存時間を、時間、時間における群における群のリスク集合の大きさをの注目しているイベント数をとする。

手順 1 群 について、時間における注目しているイベントの CIF 推定量

手順 2 群における注目しているイベントあるいは競合リスクを共にイベントとについて、時間した Kaplan-Meier 推定量
手順 3を算出する。

を算出する。
か ら、 修 正 さ れ た リ ス ク 集 合 の 大 き さ をとする。
手順 4 スコアを算出する。このとき、は漸近的に自由度 1 のカイ二乗分布に従う。
手順 4 のスコアはのように考えると log-rank 検定と類似の形式になる。と変形するとのが重みで、log-rank 検定の一般化のように他の重みを採用することもできる。手順 4の分散の導出は、Pintilie(2006, 付録 A.3) に詳しい。

数値例として、表 19 を取上げる。なお、同表は正しい結果が得られるために Pintilie(2006) の表 5.2に最後の 2 行を加えた。イベントの 0 は打切り、1 は注目しているイベント、2 は競合リスクである。
共変量は A,B の 2 群からなる。

― 70 ―

This document is provided by JAXA.

表 19 数値例

R プログラムを採用するとを得る。すなわち、より、注目するイベントの 2 群間には有意差がない。

― 71 ―

This document is provided by JAXA.

(4)CIF に関する 2 群間の差の Pepe&Mori(1993) 検定
競合リスクを伴う 2 群間の差を検定しよう。群群を合併した観測値の CIF を、個体数をとする。2について、を検定しよう。打切り、あるいは競合リスクをイベントとした生存関数の Kaplan-Meier 推定量をとしたとき、について、s は漸近的に正規分布に従う。

は自由度 1 のカイ二乗分布、すなわち、はで推定される。各群 に関する分散はとなる。ここにで、群 の時間でのリスク集合の個体数をのイベント数を、群 の時間、群 の時間の注目している、あるいは競合リスクでの注目しているイベント数をベントの CIF を F、群 の競合リスクに対する CIF を、群 における注目しているイとする。Lunn(1988) は、3 群以上の場合へ拡張している。
表 9 のマウス発癌性データについて、Peto&Mori 検定を行うとより、胸腺リンパ腫を注目するイベントとしたとき、2 群間の CIF は高度に有意になる。R プログラ
ムは Pintilie(2006, 付録 B.3.2) に公開されている。
(5) 部分分布に基づく競合リスク回帰モデル
部分分布のハザード関数 (32) 式に基づく競合リスク回帰モデル (Competing risks regression) は、共変量が 1 個の場合、で与えられる (Fine&Gray,1999;Pintilie,2006,6.2,6.3 節 )。ここには部分分布のハザード関数、― 72 ―

This document is provided by JAXA.

は部分分布のベースラインハザード関数である。このとき、尤度はとなる。ここに、は時間までに、あるイベントが起きていない個体、および時間あるいは (スクイベントが起きた個体の集合でまでに競合リかつ競合リスクイベントが起きた個体 )} となる。他のイベントが起きた個体は、すべての時間でリスク集合に残る。重みはから算出する。ただし、としたとき、イベント (の生存時間に対する Kaplan-Meier 推定量である。注目しているは打切り分布とする ) が起きた時間でリスクに入る個体 ( とする ) の集合は、 より前に競合リスクイベントが起きた ( この場合、重みは 1) 個体、および時間までにいかなるタイプのイベントも起きていない ( この場合、重みは 1 より小さい ) 個体から成る。競合リスクイベントが起きた個体は、部分尤度に加味されない。
重みの算出法は下記の通りである。

手順 1: 縦軸に注目しているイベントのみの時点 ( )、横軸にすべての個体 ( ) の番号を記入し、各セルの重みをとする。

手順 2: 各 ( ) について、となるを 1 とする。

手順 3: 注目しているイベントが起きた個体、および打切り例について、以外のセルは、部分尤度に加味しない ( × を付ける )。
手順 3: 残りのセルにはを割付ける。
(49) 式から、スコア統計量を得る。CIF の予測はから算出する。ここに、は部分分布の累積ハザードで、予測したい共変量の値をとしたとき、Breslow 型推定量を採用する。
表 15 の数値例に、共変量を加えた表 20 のデータを解析する。

は表 20 のように算出される。

― 73 ―

This document is provided by JAXA.

は、表 21 のようになる。よって、(53) 式のはと書ける。
表 20 共変量を加えたデータ

表 21 重みの計算

R プログラムから― 74 ―

This document is provided by JAXA.

が求まりより、共変量は高度に有意である。となる。
次に、表 14 のデータを解析する。MGUS を注目しているイベント、骨髄腫とその他を競合リスクとして再解析する。R プログラムを採用するとが得られる。表 22 に周辺モデルおよび層別比例ハザードモデルの結果も載せておく。

― 75 ―

This document is provided by JAXA.

表 22 競合リスク回帰モデル、周辺モデルおよび層別比例ハザードモデルとの比較

周辺モデル

共変量

競合リスク回帰(p 値 )

sex-0.3132(0.097)0.0753(<0.0001)

0.0079(0.6006)-0.0024(0.907)-0.3493(0.0253)age

0.0778(<0.0001)-0.4525(0.0176)-0.0761(0.8231)0.0043(0.993)0.0516(<0.0001)

hgb-0.00765(0.250)-0.2055(<0.0001)-0.0998(<0.3329)

0.0573(0.727)-0.1669(0.002)

mspike0.02976(0.160)0.1244(0.579)-0.6917(0.0995)-0.3071(0.603)-0.0946(0.257)

死亡多発性骨髄腫

その他

層別比例ハザード

競合リスクが存在しない場合、死因別ハザードと部分分布のハザードに基づく回帰モデルが一致することを確認しよう。表 20 のデータで、競合リスクイベントをすべて注目するイベントとした表 23 を解析する。
表 23 表 20 のデータで、競合リスクイベントをすべて注目しているイベントとしたデータ

従来の比例ハザードモデル ( 死因別ハザード ) の場合、R プログラムを採用すると― 76 ―

This document is provided by JAXA.

を得る。一方、部分分布のハザードに基づく競合リスク回帰モデルの場合、R プログラムを採用すると、を得、両者が一致することが分かる。
(6)CIF の比例ハザード性の検定
CIF の比例ハザード性の検定を行うことができる (Pintilie,2006,6.2.3 節 )。表 14 のデータについて、TIME 関数を用い、性別との交互作用を新たな共変量にすれば良い。R プログラム― 77 ―

This document is provided by JAXA.

を採用すればが得られる。すなわち、性別× time の交互作用は有意にならないため、比例ハザード性は妥当といえる。また、ハザード比は t と共にで変動する。

次に、注目しているイベントの CIF を F としたとき、の 2 次元プロッ

トで視覚的にチェックすることもできる。R プログラム

から、図 10 が得られる。

図 10 比例ハザード性の検証― 78 ―

This document is provided by JAXA.

(7) 時間依存型への拡張
時間依存型共変量を伴う場合へ競合リスク回帰モデルを拡張しよう (Beyersmann& Schumacner,2008)。例えば、図 11 のように“集中治療室 (ICU) に長く留まれば感染症肺炎のため ICU 内での死亡が増えるか”という問題を考える。患者 # 3163 は、ICU 入室中に他の死因 ( 競合リスク ) で死亡した。
# 3147 は ICU を退出したが、その後、感染症肺炎で死亡した。# 30148 は、ICU 退出後、入院中 ( すなわち、打切り ) である。# 30236 は、ICU 入室中の途中までは、異常が無かったが、ICU 退室後、しばらくして感染症肺炎を発症し、入院中である。# 30254 は、ICU 入室中、途中までは異常が無かったが、退室後しばらくして感染症肺炎で死亡した。# 1014378 は、
ICU 入室中の途中までは異常は無かったが、しばらくして感染症肺炎のため ICU 内で死亡した ( 注目しているイベント )。

図 11 ICU

表 24 のデータは、それぞれの患者の時間依存型データである。各変数はを表している。ただし、の場合は、の値を無視する。

表 24 ICU データ― 79 ―

This document is provided by JAXA.

R プログラムを採用すると、を得る。
“ICU 内感染症肺炎”発症の有無は高度に有意で、院内肺炎のハザード比は 3.016 となる。
R のパッケージ‘kmi’を確認するため、3.1 節の PBC データを解析してみよう ( この R パッケー
ジは http://cran.r-project.org/web/packages/kmi を参照 )。このデータは時間依存型であるが、注目しているイベント ( 死亡 ) のみで、競合リスクは存在しない。よって、表 25 のように書ける。
表 25 PBC データ― 80 ―

This document is provided by JAXA.

R プログラムを採用するととなり、3.1 節の結果と一致する。

参考文献
Aalen,O., Borgan, Φ , and Gjessing,H.K. (2008): Survival and Event History Analysis. Springer.
Aitkin, M., Laird, N., and B. Francis, B.(1983). A reanalysis of the Stanford heart transplant data,
Journal of the American Statistical Association 78, 264-292.
赤澤宏平，柳川尭 (2010): サバイバルデータの解析 ―生存時間とイベントヒストリーデータ―、近代科学社．
Allison,P.D.(1995):Survival Analysis Using the SAS System, A Practical Guide. SAS Publisher
Altman, D.G., & De Stavola, B.L.(1984). Practical problems in fitting a proportional hazards model
to data with updated measurements of the covariates, Statistics in Medicine, 13, 301-341.
Andelsen, P.K., and Gill, R.D.(1982). Cox’s regression model for counting process: A large sample
study, Annals of Statistics, 10, 1100-1120.
Beyersmann,J. and Schumacher,M.(2008):Time-depependent covariates in the proportionalsubdistribution hazards model for competing risks, Biostatistics, 9, 765-776.
Collet,D.(2003):Modelling Survival Data in Medical Research, 2nd edi. Chapman & Hall.
Cox,D.R.(1975): Partial likelihood, Biometrika, 62,269-276.
Christensen, E., Schlichting, P., Andersen, P.K., Fauerholdt, L., Schou, G., Pedersen, B.V., Juhl, E.,
Poulsen, H., and Tygstrup, N., Copenhagen Study Group for Liver Disease (1986). Updating
prognosis and therapeutic effect evaluation in cirrhosis with Cox’s multiple regression
model for time-dependent variables, Scandinavian Journal of Gastroenterology, 21, 163-174.
Fine,J.P., and Gray,R.J.(1999): A proportional hazards model for the subdistribution of a compet― 81 ―This document is provided by JAXA.

ing risk, Journal of the American Statistical Association, 94, 496-509.
Gray, R.J.(1988): A class of k-sample tests for comparing the cumulative incidence of a competingrisk, Annals of Statistics, 16, 1141-1154.
Gray,R.J.(2010):cmprsk:Subdistribution analysis of competing risks. R package version 2.2-1,
URL http://cran.r-project.org/web/packages/cmprsk.
Kalbﬂeisch,J.D.,and Prentice,R.L.(2002):The Statistical Analysis of Failure Time Data. John Wiley.
Klein,J.P.,and Moeschberges,M.L.(2003): Survival Analysis, Springer.
Lee,E.M., et al.,(2003):Statistical Methods for Survival Data Analysis. John Wiley.
Lunn,M.(1998): Applying k-sample tests to conditional probabilities for competing risks in aclinical trial, Biometrics, 54, 1662-1672.
中村剛 (2001):Cox 比例ハザードモデル , 朝倉書店
西川正子 (2008): 生存時間解析における競合リスクモデル , 計量生物学 , 29,141-170.
大橋靖雄 , 浜田知久馬 (1995): 生存時間解析 , 東大出版
Pepe, M.S., and Mori,M.(1993): Kaplan-Meier, Marginal or conditional probability curves in summarizing competing risks failure time data, Statistics in Medicine, 12, 737-751.
Pintilie,M.(2006):Copeting Risks, a Practical Perspective. John Wiley.
Prentice,R.L., Williams,B.J., and Peterson, A.V.(1981): On the regression analysis of multivariate
failure time data, Biometrika, 68, 373-379.
Putter,H., Fiacco, M., and Geskus, R.B.(2007): Tutorial in Biostatistics; Competing risks andmulti-state models, Statistics in Medicine, 26, 2389-2430.
Therneau,T., and Grambsch,P.(2000): Modeling Survival Data: Extending the Cox Model.
Springer-Verlag.
Tsujitani, M. and Koshimizu,T(2000): Neural Discriminant Analysis, IEEE Transaction on NeuralNetworks, 11, 1394-1401.
辻谷将明 , 外山信夫 (2007):R によろ GAM 入門 , 行動計量学 ,34,111-131.
辻谷将明 , 竹澤邦夫 (2009):R で学ぶデータサイエンス 6 マシンラーニング , 共立出版 .
辻谷将明 , 和田武夫 (2012):R で学ぶ確率・統計 , 共立出版 .
Tsujitani,M. and Sakon,M.(2009): Analysis of survival data having time-dependent covariates,
IEEE Transaction on Neural Networks, 20,389-394.
Tsujitani, M., and Tanaka, Y.(2011):Cross-validation,Bootstrap, and support vector machines, Advances in Articial Neural Systems, 2011, Article ID302572, 6pages.
Tsujitani, M., and Baesens, B.(2012) : Survival analysis for personal loan data using generalized
additive models, Behaviormetrika, 39,9-23.
Tsujitani,M., Tanaka,Y., and Sakon,M.(2012):Survival analysis with time-dependent covariates
using generalized additive models, Computational and Mathematical Methods in Medicine, 2012,
Article ID986176, 9 pages.
Tsujitani, M.,Iba,K., and Tanaka,Y.(2012):Neural discriminant models, bootstrapping andsimulation, Computational and Mathematical Methods in Medicine, 2012, Article ID820364, 12pages.
Tsujitani, M., and Tanaka,Y.(2013):Analysis of heart transplant survival data using generalized
additive models, Computational and Mathematical Methods in Medicine, 2013, Article ID609857,
7 pages.
― 82 ―

This document is provided by JAXA.

Wood, S. N. (2000): Modelling and smoothing parameter estimation with multiple quadratic
penalties. Journal of the Royal Statistical Society, Series B, 62, 413-428.
Wood, S. N. (2006): Generalized additive models: An Introduction with R. Chapman & Hall / CRC.
Wood, S. N. (2008): Fast stable direct fitting and smoothness selection for generalized additivemodels, Journal of the Royal Statistical Society Series B, 70, 495-518.

― 83 ―

This document is provided by JAXA.

This document is provided by JAXA.