統計数理（2020）
第 68 巻 第 1 号 147–174©2020 統計数理研究所

特集「複雑な依存構造を持つデータの統計解析法 — コピュラとその周辺 — 」［総合報告］

コピュラを用いた生存時間解析——相関のあるエンドポイントとメタ分析の活用——

江村 剛志1 ・道前 洋史2（受付 2019 年 4 月 24 日；改訂 2020 年 1 月 18 日；採択 1 月 20 日）

要旨

ここ 10 年間で医学・生物学関連データベースの公開が急速に発展し，個々の患者の詳細かつ正確な情報が誰にでも利用出来るようになった．これらデータベースでは，患者の生存期間，無増悪期間，腫瘍径，遺伝子発現量など患者の極めて詳細な情報が記録されている．加えて，複数のデータソースを統合して解析する「メタ分析」の重要性がここ 10 年間で高まり，複雑かつ膨大な生存期間データを古典的な生存時間解析法だけで解析することは不十分になりつつある．本稿では，コピュラを用いて 2 つの生存期間変数をモデリングする手法について総説する．モデルの実践的な有用性を示すため，より具体的に患者の生存期間と無増悪期間の相関をコピュラでモデル化する手法を考える．また複数のデータソースを統合して，患者の生存期間と術後無増悪期間を同時解析する「メタ分析」
法のためのモデルを紹介する．モデルのパラメータを推定する際に，生存期間と無増悪期間が半競合リスク関係にあることを考慮した上で最尤法を用いる必要性を解説する．最後に，個別医療の問題に関連した応用例の一つとして，遺伝子発現量と増悪情報を用いて死亡率の動的予測を行う手法を紹介する．
キーワード：遺伝子発現量，個別医療，競合リスク，動的予測，臨床試験，Cox 比例
ハザードモデル．

1.

はじめには，最も歴史のある学術分野の 1 つでその起源は生命表が生存時間解析
（Survival analysis）
考案された 1600 年代にまで遡る．生存時間解析で使用されている統計手法は極めて多岐に及び，いまなお多くの研究者が学術研究を行っていることは生存時間解析のハンドブック（Kleinからも理解できる．生存時間解析は現実問題にも幅広く応用されており，医学研究et al., 2014）における癌患者の生存率の計算法，工場で生産される部品の信頼性の計算法，政府統計または保険数理における生命表の計算法など多くの統計手法の数理的基礎を成す学術分野といえる．
医学研究における生存時間解析では，年齢や性別などの共変量が生存期間にどのように影響で調べることが頻繁に行われている．Cox 回帰法は生存時するのかを Cox 回帰法（Cox, 1972）
間解析で最も広く使用されてきた統計手法の一つであり，これからも医学研究に必要不可欠な手法であることは疑う余地が無い．
（癌ゲノムのデータベー
しかしながら，ここ 10 年間で医学・生物学関連の公開データベース
1

長庚大学 情報管理学系：〒 333 桃園市亀山区文化一路 259 号（台湾）；takeshiemura@gmail.com

2

北里大学 薬学部：〒 108–8641 東京都港区白金 5–9–1

148

統計数理

第 68 巻 第 1 号 2020

図 1．術後 t 日目から開始する患者の死亡率の動的予測．

ス GEO や TCGA など）
や統計解析フリーソフト（R など）が急速に発展し，個々の患者の詳細かつ正確な情報が誰にでも利用出来るようになった．これらデータベースでは，患者の生存期間，無増悪期間，腫瘍径，遺伝子発現量など患者の極めて詳細な情報が記録されている．加えて，複数のデータソースを統合して解析する「メタ分析」の重要性がここ 10 年間で高まり，複雑かつ膨大な生存期間データを古典的な生存時間解析法だけで解析することは不十分になりつつある．このような患者の個別データを有効利用し，個別医療の進展に貢献できる統計手法を開発することは生物・医学統計学者の重要な使命でもある．
本稿では，コピュラを用いて 2 つの生存期間変数をモデリングする手法について総説する．
モデルの実践的な有用性を示すため，より具体的に患者の生存期間と無増悪期間の相関をコ
ピュラでモデル化する手法を考える．また複数のデータソースを統合して，患者の生存期間と無増悪期間を同時解析する「メタ分析法」のためのモデルを紹介する．モデルのパラメータを推定する際，生存期間と無増悪期間が半競合リスク関係（Fine et al., 2001）にあることを考慮した上で最尤法を用いる必要性を解説する．最後に，個別医療の問題に関連した応用例の一つとして，遺伝子発現量と増悪情報を用いて死亡率の動的予測を行う手法を紹介する（Emura et al.,2018）．すなわち，患者個々の術後経過と遺伝子情報を考慮した死亡率の個別予測式をメタ分
析のモデルから導出し，患者レベルの死亡予測を行う手法（図 1）を紹介する．このような動的予測の手法は多く提案されている（Sène et al., 2016; Proust-Lima et al., 2014; Rondeau et al.,が，本稿で紹介する手法はメタ分析を使用する点と，コピュラを利用する点で既出2017 など）のものと大きく異なる．
本稿は以下のように構成されている．2 節では生存時間解析における重要な用語「エンドポイント」を解説し，なぜコピュラを利用するのかを癌の臨床研究の側面から説明する．3 節では，生存時間解析におけるコピュラのモデルを定式化し，提案されているいくつかのモデルと統計手法を紹介する．4 節では Frailty モデルを用いたメタ分析の考え方を紹介する．5 節ではFrailty モデルにコピュラを組み合わせた統計モデルを紹介し，エンドポイント間の半競合リスク関係を解説する．6 節では個別医療の問題に関連した動的予測の手法を紹介する．5 節，6 節
では卵巣癌患者のデータを用いた事例を用いて，統計手法の利用法を解説する．
2.

生存期間のエンドポイント

医学研究における生存時間解析の基本的な目的の 1 つは，特定の治療や予後因子が，患者の生存期間にどれだけ影響するか調べることである．医学研究において，生存期間はエンドポイ
ントと呼ばれ，具体的には次のような複数の定義がある．
• 全生存期間（OS: Overall survival）

コピュラを用いた生存時間解析

149• 無増悪期間（TTP: Time to tumor progression）• 無増悪生存期間（PFS: Progression-free survival）• 無病生存期間（DFS: Disease-free survival）

増悪とは癌の再発，腫瘍増大，遠隔転移などを意味する．無増悪期間とは増悪がなかった期間である．例えば，手術後に残った直径 1cm の腫瘍がそのまま増大せず，かつ遠隔転移が見られなければ，無増悪とみなされる.
エンドポイントの定義は上記以外にもあるが，いずれにせよ定義が明確であり，かつ臨床的に重要な転帰を表すものでなければならない（Pazdur, 2008; Piedbois and Croswell, 2008; Le．医学研究者や臨床医は，これらエンドポイントから重要なものを選び，Tourneau et al., 2009）
個々の患者に対して追跡・観察する．
治療効果や予後因子がこれらエンドポイントにどのように影響するのかを調べるためには通を用いる．例えば，生存期間を応答変数，治療群と予常 Cox 比例ハザードモデル（Cox, 1972）
後因子を共変量とし，部分尤度法により共変量がハザード（瞬間死亡率）に与える影響を調べることが出来る．このような共変量の効果を調べる解析は，医学研究において極めて頻繁に行われている．
全生存期間は癌の研究において最も重要かつ明確に定義できるエンドポイントであり，真またはゴールドスタンダードエンドポイント（Pazdur, 2008;のエンドポイント（true endpoint）と呼ばれる．例えば，頭頸部癌のShi and Sargent, 2009; Michiels et al., 2009; Oba et al., 2013）
化学療法や放射線療法の効果を評価する際の最も重要なエンドポイントは全生存期間である（Michiels et al., 2009; Le Tourneau et al., 2009）．
無増悪期間，無増悪生存期間，無病生存期間などのエンドポイントはしばしば代替エンドポとよばれ，全生存期間に次いで重要なエンドポイントである．代替イント（surrogate endpoint）
エンドポイントを使用する場合，“増悪” の判定は臨床医の恣意性が伴うため，RECIST ガイドなどに従ってすすめることが推奨されている．
ラン（Eisenhauer et al., 2009）（例えば，全生存期間と無増悪生存期
多くの臨床研究では，患者から 2 つのエンドポイント間）を観察するためのフォローアップを行う．その後，得られたデータをもとに，Cox 回帰法などで治療効果や予後因子が各エンドポイントに与える影響を解析する．しかしながら，単純な Cox 回帰では 1 つのエンドポイントを応答変数とした周辺解析しか行うことができず，2 つ以上のエンドポイントを同時解析する場合には不十分である．
例えば，ある治療が生存期間と無増悪期間の両方に与える影響に興味があるとする．この場合，生存期間を応答変数とした Cox 回帰と，無増悪期間を応答変数とした周辺 Cox 回帰を別々に実行することが考えられるが，これらの 2 つの解析結果からは，全生存期間と無増悪期間の相関の推定値は得られない．無増悪期間が全生存期間とどのように関連しているかを調べ
るには相関の解析が必要であり，このような相関の情報は抗癌剤開発の観点からも重要である．加えて，無増悪期間が全生存（Sherrill et al., 2008; Michiels et al., 2009; Rondeau et al., 2015）
期間の代替エンドポイントであることを実証したい場合や，無増悪期間を全生存期間の予後予測に利用する場合に，相関の推定値が必要であることを考慮すると，周辺 Cox 回帰法はデータを有効に活用しているとはいえない．
3.

コピュラの生存時間モデル

コピュラは 2 つの確率変数間の相関構造をあらわす関数である．コピュラの概念は数学者
Abe Sklar により導入された（Sklar, 1959）．コピュラの数学的な定義と性質は Nelsen（2006）の著書に良くまとめられている．これは数学書であるが，大学初年度の微積分の知識で読める良

150

統計数理

第 68 巻 第 1 号 2020書である．以降では，コピュラを生存時間モデルに応用するための基礎知識や用語を概説する．
3.1 2 次元生存時間モデル
コピュラは 2 つの生存期間変数間の相関構造をモデリングするために用いることが出来る．いま X と Y を生存期間変数とし，Z1 = (Z11 , . . . , Z1p1 ) と Z2 = (Z21 , . . . , Z2p2 ) を X とY それぞれに関連する共変量とする．ここで Z1 と Z2 は Z1 = Z2 であっても良い．また
SX (x|Z1 ) = Pr(X > x|Z1 ) と SY (y|Z2 ) = Pr(Y > y|Z2 ) を X と Y の生存関数とする．このとき，X と Y の 2 次元生存関数を次のように定義する．
(3.1)

Pr(X > x, Y > y|Z1 , Z2 ) = Cθ {SX (x|Z1 ), SY (y|Z2 )}.

ここで Cθ はコピュラ（Sklar, 1959; Nelsen, 2006）で，θ は X と Y の相関の強さを定めるパラメータである．コピュラの種類は無数にあるが，本稿で主に扱うコピュラは次の 3 つである．
• 独立型: C(u, v) = uv, 0 ≤ u, v ≤ 1,
• クレイトン型（Clayton, 1978）:Cθ (u, v) = (u−θ + v −θ − 1)−1/θ ,θ > 0,

0 ≤ u, v ≤ 1,

• グンベル型（Gumbel, 1960）:1Cθ (u, v) = exp[−{(− log u)θ+1 + (− log v)θ+1 } θ+1 ],θ ≥ 0,

0 ≤ u, v ≤ 1.

独立型のコピュラを選ぶと，Pr(X > x, Y > y|Z1 , Z2 ) = SX (x|Z1 )SY (y|Z2 ) となり，X と Yは共変量を所与したときの条件付独立なモデルとなる．コピュラの定義により，上記 3 つの Cθ
はどれもコピュラであるから，Cθ は 2 次元分布関数であり，その周辺分布は [0, 1] 上の一様分布である．実際，周辺分布は Cθ (u, 1) = u と Cθ (1, v) = v であり，これらは [0, 1] 上の一様分布の分布関数である．
Cθ が 2 次元分布関数であるから，Pr(U ≤ u, V ≤ v) = Cθ (u, v) となる一様分布の確率変数の−1(U |Z1 )ペア (U, V ) を考えることが出来る．このとき，確率変数のペア (X, Y ) を変換 X = SX
−1と Y = SY (V |Z2 ) で定義すると，(X, Y ) の分布は式（3.1）を満足することがわかる．これは，確率変数のペア (U, V ) が相関構造のみを定め, SX (·|Z1 ) と SY (·|Z2 ) が周辺構造のみを定めていることを意味する．この周辺構造と相関構造が分離できることは統計モデリングの観点から極めて有用である．
（Archimedean
上記の 3 つのコピュラの例は，どれも下記で定義するアルキメデス型コピュラに属することがわかる．いま，関数 φθ : [0, 1] → [0, ∞] を連続な減少関数で φθ (1) = 0
copula）
を満足するとする．このとき，アルキメデス型コピュラはCθ (u, v) = φ−1
θ {φθ (u) + φθ (v)},で定義される．ここで φθ は生成素
（generator）と呼ばれ，その逆関数を φ−1で表す．クレイθ−θ
トン型は φθ (t) = (t − 1)/θ，θ > 0 である．その極限は limθ→0 φθ (t) = − log(t) であり，これは独立型の生成素になっていることがわかる．グンベル型は φθ (t) = {− log(t)}θ+1 で生成され，θ = 0 が独立型の生成素になる．すなわち，クレイトン型とグンベル型では極限がlimθ→0 Cθ (u, v) = uv となっているので，独立型に帰着する．表 1 にコピュラとその生成素をまとめた．
3.2 ケンドール順位相関係数
ケンドール順位相関係数（以下，ケンドール相関と記す）
は X と Y の相関を表す指標である．

コピュラを用いた生存時間解析

151

表 1．コピュラの例．

図 2．500 組の (U, V ) をクレイトン型コピュラから生成し，散布図にしたもの．

モデル（3.1）のもとで，X と Y のケンドール相関はコピュラ Cθ のみに依存することがわかり，

 1 1Cθ (u, v)dCθ (u, v) − 1

τθ = 400と書ける．ケンドール相関は周辺分布の構造 SX (x|Z1 ) と SY (y|Z2 ) に依存しない．よってコ
ピュラ Cθ の型を定めれば，τθ はパラメータ θ によって定まる．アルキメデス型コピュラでは，
ケンドール相関は生成素を使って次のように計算できる．

 1

τθ = 1 + 4
0φθ (t)dt,φ̇θ (t)φ̇θ (t) = dφθ (t)/dt.

この式を用いると，クレイトン型では τθ = θ/(θ + 2)，グンベル型では τθ = θ/(θ + 1) となり，ともに無相関 τ0 = 0 から正の完全相関 τ∞ = 1 までとりうる．
図 2 は 500 組の (U, V ) をクレイトン型から生成し，散布図にしたものである．図から U と Vとの間に正の相関があることがわかる．相関の程度は θ = 2（τθ = 0.5）よりも θ = 8（τθ = 0.8）が強いことがわかる．また，U ≈ 0 と V ≈ 0 付近でより強い相関を示し，これは下側裾従属と呼ばれるクレイトン型コピュラが有する特徴的な現象である．一性（lower tail dependence）
方，グンベル型では U ≈ 1 と V ≈ 1 付近でより強い相関を示し，上側裾従属性（upper taildependence）
が現れる．
−1(U |Z1 ) と Y = SY−1 (V |Z2 ) で変換すモデル（3.1）のもとで X と Y を生成するには，X = SX
−1
−1
れば良い．変換 SX と SY は減少関数であるため，変換により裾従属の位置は反転する．すなわち，クレイトン型では X と Y 間に上側裾従属，グンベル型では X と Y 間に下側裾従属があらわれる．

統計数理

152

第 68 巻 第 1 号 2020

3.3 オッズ比による相関
生存時間解析においては，2×2 クロス集計表と関連した相関の指標が，推定や予測の局面で有用となることがある．例えば，死亡と増悪の関連をみるとき，死亡と増悪のあり・なしを
2×2 クロス集計表に集計し，そのオッズ比を計算して相関の度合いを調べることがある．このが提案した 2 次元生命表に遡り，Oakes（1989）がより明確
オッズ比の考え方は Clayton（1978）に定式化した．
（3.1）を仮定する．オッズ比関数（英語
いま X を無増悪期間 ，Y を全生存期間とし，モデルを次のように定義する．
では Cross-ratio function）
Pr(X = x, Y = y|Z) Pr(X > x, Y > y|Z)
= Rθ {SX (x|Z1 ), SY (y|Z2 )}Pr(X = x, Y > y|Z) Pr(X > x, Y = y|Z)

左辺はオッズ比の定義であり，表 2 の 2×2 クロス集計表に基づく．右辺は，モデル（3.1）から導出された結果であり，[1,1]Cθ (u, v)Cθ (u, v)
Rθ (u, v) = [1,0][0,1]Cθ (u, v)Cθ (u, v)
はコピュラのオッズ比関数である．ここで
[1,0]

Cθ

(u, v) =∂Cθ (u, v),
∂u

[0,1]

Cθ

(u, v) =∂Cθ (u, v),∂v

はコピュラの条件付分布関数である．すなわち U = u の条件下の V の分布関数は Pr(V ≤
[1,0]
[0,1]
v|U = u) = Cθ (u, v) であり，V = v の条件下の U の分布関数は Pr(U ≤ u|V = v) = Cθ (u, v)である．また，∂2[1,1]
Cθ (u, v) =Cθ (u, v)∂u∂v
はコピュラの密度関数を表す．
関数 Rθ (u, v) は点 (u, v) における局所的な相関を表し，次のように解釈できる．
• Rθ (u, v) > 1；局所正相関，• 0 < Rθ (u, v) < 1；局所負相関，• Rθ (u, v) = 1；局所独立．

独立型コピュラのもとで Rθ (u, v) = 1 が全ての 0 ≤ u ≤ 1 と 0 ≤ v ≤ 1 で成立する．クレイト
ン型コピュラではオッズ比関数が定数 Rθ (u, v) = 1 + θ となり，これは原論文 Clayton（1978）に由来する．アルキメデス型コピュラのオッズ比関数は次のように書けるRθ (u, v) = rθ {Cθ (u, v)},rθ (s) = −sφ̈θ (s)/φ̇θ (s),φ̈θ (t) = d2 φθ (t)/dt2 .

オッズ比関数 Rθ が相対リスク（Relative risk）
として解釈できることは生存時間解析において極めて有用である．例えば，X を無増悪期間 ，Y を全生存期間とし，増悪が生存期間に与える影響に興味があるとする．ある時刻 x において増悪したか否かの情報を共変量と考え，次の 2
つの条件付ハザード関数を定義する．
表 2．2 × 2 クロス集計表．A, B, C, D はイベントが起こる頻度，オッズ比は AD/(BC)．

コピュラを用いた生存時間解析

153• λY (y|X = x, Z) = Pr(y ≤ Y < y + dy|Y ≥ y, X = x, Z)/dy

時刻 x で増悪したときの条件付ハザード，• λY (y|X > x, Z) = Pr(y ≤ Y < y + dy|Y ≥ y, X > x, Z)/dy
時刻 x 以前に増悪しなかったときの条件付ハザード．
相対リスクは 2 つの条件付ハザード関数の比で定義される；
λY (y|X = x, Z)= Rθ {SX (x|Z1 ), SY (y|Z2 )}.
λY (y|X > x, Z)

右辺を導出するには，モデル（3.1）の仮定が必要である．Rθ > 1 であれば，増悪したことが死亡リスクを上昇させることを意味する．クレイトン型コピュラでは，相対リスクが定数 1 + θ なので，これは比例ハザード性を意味する．また，θ > 0 なので，増悪と死亡が任意の時刻 (x, y)
において局所正相関をもつモデルを表す．
Clayton（1978）
は，データを表 2 の 2×2 クロス集計表にまとめ，そのオッズ比を計算することにより θ を推定することを提唱している．この手法は X と Y が共に右側打ち切りを受けている場合にも応用でき，さらに周辺分布の構造 SX (x|Z1 ) と SY (y|Z2 ) に依存しないセミ・パラ
メトリックな推定法である．この手法はアルキメデス型コピュラに一般化可能であり（Emuraet al., 2010），さらに半競合リスクデータ（Wang, 2003）や従属切断データ（Emura and Wang,にも適用可能である．
2010; Emura et al., 2011）
3.4 応用例
コピュラは生存時間のモデルに広く用いられてきたが，その中でもとりわけ臨床医学研（two-step
究において広く利用されているのが Burzykowski et al.（2001）による 2 段階推定法
である．代替エンドポイントを用いる臨床試験の策定にあたり，現時点ではコピュラ
method）
を用いた 2 段階推定法（Burzykowski et al., 2005）が広く用いられている．最近，この 2 段階推定法を実行する R パッケージ（Rotolo et al., 2018）が開発され，利用が容易になった．この Rパッケージではクレイトン型，グンベル型，プラケット型コピュラが利用できる．
いま，i 番目の研究単位から得られた j 番目の患者に対して，Xij を無増悪期間，Yij を全生
は Xij と Yij の周辺ハザードに存期間とする．Burzykowski et al.（2001）



λXij (t|Zij ) = λXi (t) exp(βXi Zij )λY ij (t|Zij ) = λY i (t) exp(βY i Zij )のような Cox モデルを仮定した．ここで Zij は治療群を表す変数，β· はその回帰係数，λ·i (·)
は i 番目の研究単位の基準ハザード関数である．さらに，Xij と Yij の同時分布を，Pr(Xij > x, Yij > y ) = Cθ

  x

exp

−



λXij (t|Zij )dt
0

  y

,

exp

−



λY ij (t|Zij )dt0で定義し，Cθ はコピュラとした．パラメータ推定の方法はいくつかある．例えば，まず通常の Cox 回帰法（部分尤度法）やワイブル回帰法で (β· , λ·i (·)) を推定してから，次に θ を推定する手法がある．他にも，(θ, βXi , βY i , λXi (·), λY i (·)) を同時推定する全尤度法がある．さらに基準
ハザードの共通性 λXi (·) = λX0 (·) と λY i (·) = λY 0 (·) を仮定して推定する場合もある．いずれにせよ，独立な右側打ち切りのある生存期間データに適用可能である．詳細は Burzykowski et，Rotolo et al.（2018）
を参照されたい．
al.（2001, 2005）
上記手法において，エンドポイント間の相関を考慮しているにもかかわらず推定の段階において，独立右側打ち切りしか想定していない点に注意すべきである．例えば，無増悪期間を

統計数理

154

第 68 巻 第 1 号 2020

Xij として βXi を Cox 回帰法（部分尤度法）で推定する場合，死亡を独立右側打ち切りとみなして解析していることになる．しかし死亡による打ち切りは，無増悪期間と明らかに独立でないため，多くの患者が早期に（増悪する前に）死亡している場合には回帰推定値に大きなバイアスが入る．この問題は従属打ち切りの問題として知られ，打ち切り数が多いほど，また相関が強いほどバイアスが大きくなることが知られている（Emura and Chen, 2016, 2018）．
この従属打ち切り問題を統計的に正しく扱うには，データ構造を半競合リスク（Fine et al.,の枠組みで捉えた統計手法を用いると良い．すなわち，死亡は増悪の観測を打ち切るた
2001）
め競合リスクであるが，増悪は死亡の観測を打ち切ることがないため競合リスクでないとする考えである．半競合リスクデータの枠組みで正しく回帰推定値を得る場合には，部分尤度法（Cox 回帰法）ではなく，全尤度を利用した推定法（Chen, 2012; Emura et al., 2017）を利用することになる．半競合リスクデータを扱う統計手法は 5 節で再考する．
上記以外にも，コピュラは極めて多くの生存時間解析法で利用されてきた．例えば，クラスター生存時間解析（Rotolo et al., 2013; Emura et al., 2017; Peng et al., 2018），競合リスク
解析（Lo and Wilke, 2010; de Uña-Álvarez and Veraverbeke, 2013, 2017; Shih and Emura, 2018;
Shih et al., 2019; Zhou et al., 2019; Emura et al., 2019c; Wang et al., 2020），半競合リスク解析（Wang, 2003; Chen, 2012; Peng, 2019; Peng and Xiang, 2019），従属打ち切り下での生存時
間解析（Emura and Michimae, 2017; Moradian et al., 2019; Emura and Chen, 2016, 2018），従属切断下での生存時間解析（Chaieb et al., 2006; Emura and Wang, 2012; Emura and Murotani,
2015; Emura and Pan, 2020），ヴァインコピュラを用いた多変量生存時間データ（Barthel et al.,
，繰り返しイベント時間解析などである（Ling et al., 2016; Li et al., 2019, 2020）．
2018）
4.

メタ分析の活用

メタ分析とは，複数の研究単位から集められたデータを統合し，より一般的な母集団に対する推論を行う手法である．メタ分析は出版された複数の研究論文から要約統計量（相対リスクなど）を抽出して統合解析を行う手法を指すこともあるが，本稿でいうメタ分析とは患者個人（Individual-patient data meta-analysis）を指す．複のデータを用いる，いわゆる IPD メタ分析
数の研究単位から集めた患者の個別データに基づく統合解析は，全生存期間と代替エンドポイ
ントの相関解析に不可欠とされている（Burzykowski et al., 2001, 2005; Shi and Sargent, 2009;．
Rotolo et al., 2018; Sofeu et al., 2019, 2020）
IPD メタ分析を通じて全生存期間とその代替エンドポイントの相関を推定する研究が数多，頭頸部癌（Michiels et al., 2009），胃癌（Oba etく行われている．大腸癌（Buyse et al., 2008），結腸癌（Alonso and Molenberghs, 2008）の患者において全生存期間と無増悪生存期
al., 2013）
間の間に強い正相関があることが報告されている一方，乳癌の患者においては正相関がそれほ．このように全生存期間と無増悪生存ど強くないことが報告されている（Michiels et al., 2016）
期間の間に意味のある相関があるかどうかは癌種や進行度などにより異なる．これらのメタ分析において，相関の推定にはコピュラが用いられてきた（3.4 節）．
推定効率の点からも，2 つ以上のエンドポイントをコピュラでモデル化した同時解析には利は，連続変量の多変量メタ分析において，同時モデルで推定点があろう．実際，Riley（2009）
した回帰係数が，周辺モデルで推定したものより推定精度が良く（推定量の平均二乗誤差が小さく）なるとしている．
ほとんどのメタ分析においては，研究単位間の異質性を考慮した解析が行われる．一般的に，異なる時期，異なる地域・国で行われた研究を統合する場合，試験間で，参加者の背景，
エンドポイントや予後因子の定義など，さまざまな要因が厳密には異なることが一般的であコピュラを用いた生存時間解析

155

図 3．卵巣癌患者を 4 つのデータソースから集めたメタ分析のデータ（Ganzfried et al.,2013; Emura et al., 2018）．データは R パッケージ joint.Cox（Emura, 2019）で入手可能．

る．また，設備が整い経験豊富なスタッフが揃った国は，そうでない国より，患者の死亡や増悪リスクが低いと考えられる．ほとんどのメタ分析の例において，このような研究単位間の異質性は，患者から観測されている共変量で説明しきれないことが多く，これを観測不能なランダム効果として解析モデル内に取り入れる．このような解析モデルをランダム効果モデルと呼ぶ．
（Random-eﬀects model）
図 3 は 4 つの研究単位から集められた患者をもとにしたメタ分析のためのデータの例を示している．データは合計 912 人の外科手術後の卵巣癌患者を追跡調査し，増悪があったかどうかを調べたものである．図 3 から増悪率が研究単位間で大きく異なることがわかる．例えば，最，最も低いもので 49% である（Study 4）．この増悪率の差異はも高い研究単位で 83%（Study 2）
追跡期間にもよるため，Study 2 が Study 1 より増悪率が高いとは断言できず，増悪率の異質性はデータから容易に観測できるものではない.
生存時間解析において，ランダム効果モデルに相当するものが，Frailty モデルである．とり（Share
わけ Shared frailty といわれるモデルにおいては，研究単位内の全ての患者で共通したされた）観測されない因子が存在することを仮定し，この因子の分散パラメータで異質性の度合いをモデル化する（Duchateau et al., 2002; Rondeau et al., 2015）．
図 3 のデータを解析する統計手法を簡単に紹介する．G 個の研究単位を考え，i 番目の研究単位が Ni 人の患者を含むとする（i = 1, 2, . . . , G）．i 番目の研究単位の j 番目の患者に対し，Xij をイベント時間，Zij を共変量ベクトルと定義する．図 3 の例で（ j = 1, 2, . . . , Ni ）は，G = 4 で N1 = 84, N2 = 58, N3 = 260, N4 = 510 であり，合計 912 人となる．
Shared frailty モデルでは，i 番目の研究単位の j 番目の患者（ j = 1, 2, . . . , Ni ）のハザード関
数を次のようにモデル化する
λij (t|ui , Zij ) = ui λ0 (t) exp(β  Z ij ).

こ こ で β = (β1 , . . . , βp ) は 未 知 の 回 帰 係 数 ，λ0 (·) は 未 知 の 基 準 ハ ザ ー ド 関 数 ，ui > 0（i = 1, 2, . . . , G）は観測不能な Frailty 項である．上記モデルにおいて，Frailty 項 ui が全ての患者で共有されていることが Shared frailty と呼ばれる由来である．
Frailty 項 ui > 0 は確率変数の実現値とみなす．以下では，最も広く使用される平均，分散のガンマ分布を用いる．すなわち，確率密度関数が



(4.1)

1
u
fη (u) =u1/η−1 exp −ηΓ(1/η)η 1/η,

η > 0,u>0であると仮定する．他の研究単位と比較して高リスクな研究単位は ui > 1，低リスクな研究単

156

統計数理

第 68 巻 第 1 号 2020

位は 0 < ui < 1 で与えられ，その分散 η は研究単位間の異質性の度合いを決めるパラメータである．極限は η → 0 研究単位間に異質性が無いことを意味する．実現値 ui は観測することが出来ないが，その分散 η はデータから推定することが出来る．
パラメータ (η, β, λ0 (·)) を最尤法で推定する様々な方法が提案されており，そのうち Frailty
分布としてガンマ分布や対数正規分布を使用する場合は，R のパッケージ等を利用して最尤推定量の計算が可能である（Hirsch and Wienke, 2012; Rondeau and Gonzalez, 2005）．
実データではイベント時間 Xij は右側打ち切りにより観測出来ないケースが多い．そのため，右側打ち切り時間 Cij を定義し，観測される時間を Tij = min(Xij , Cij ) とする．そのため
i 番目の研究単位の j 番目の患者に関するデータは (Tij , δij , Zij ) となり，δij = I(Tij = Xij ) は
イベントが起こったか否かを示す指標である．
パラメータ (η, β, λ0 (·)) の推定はデータ {(Tij , δij , Zij ); i = 1, 2, . . . , G; j = 1, 2, . . . , Ni } が与えられたときの尤度関数を構成し，対数尤度関数を (η, β, λ0 (·)) に関して最大化することによって行う（詳しくは Emura et al., 2019b の 2.5 節を参照）．
5.

メタ分析による 2 つのエンドポイントの同時解析

本節では 2 つのエンドポイントを同時に解析するための方法を紹介する．とりわけ，無増悪期間と全生存期間を同時解析し，その関連を調べるために患者の個別データ（Individualを用いたメタ分析を考える．このような分析法は「IPD メタ分析」と呼ばれPatient-Data, IPD）ており，要約データに基づくメタ分析（個別データを用いない分析）と本質的に異なる．3 節で示したようにエンドポイント間の相関をコピュラでモデル化し，4 節で示したように研究単位間の異質性を Frailty 項でモデル化するために Joint frailty-copula モデル（Emura et al., 2017）を導入する．
5.1 エンドポイントの測定と半競合リスク
複数のエンドポイントを扱う統計手法では，エンドポイントがどのような打ち切りを受けているか明確にする必要がある．これは，データ構造を明らかにするだけでなく，尤度関数を正しく構成するために必要なプロセスである．
G 個の研究単位を考え，i 番目の研究単位が Ni 人の患者を含むとする（i = 1, 2, . . . , G）．i 番（ j = 1, 2, . . . , Ni ）に対し，目の研究単位の j 番目の患者• Xij = 無増悪期間，• Dij = 全生存期間，• Cij = 右側打ち切り期間（独立かつ無情報な打ち切りを仮定），と定義する．患者は 3 つの変数 Xij ，Dij ，Cij を有しているが，どれが観測されるかは，患者の状態や打ち切り期間によって決まる．例えば，患者の状態が安定していれば，Xij と Dij は観測されず，打ち切り期間 Cij のみが観測される．この場合，打ち切り期間 Cij までに増悪と死亡が起こらなかった情報が記録される．
を示す．ケース A では，患者が増悪し，その後
図 4 で患者がとりうる 4 つのケース（A〜D）
死亡しているので，Xij と Dij が観測される．ケース B では，患者は増悪したが，死亡する前に打ち切りとなっているので，Xij と Cij が観測される．すなわち，Dij が Cij によって右側打ち切りされている．ケース C では，患者が増悪せずに死亡しているので，Dij のみが観測される．この場合，Xij が Dij によって右側打ち切りされている．Xij と Dij に相関があるので，
この打ち切りは通常の「独立打ち切り」でないことに注意する．ケース D では，Xij も Dij も観コピュラを用いた生存時間解析

157

図 4．患者がとりうる 4 つのケース（A〜D）．観測されるイベントは●，観測されないイベン
トは○で表す．観測される期間は実線，観測されない期間は点線で表す．

測されず，打ち切り期間 Cij のみが観測される．すなわち，Xij も Dij も右側打ち切りされている．
まとめると i 番目の研究単位の j 番目の患者（ j = 1, 2, . . . , Ni ）に対して観測されるデータは∗(Tij , Tij∗ , δij , δij, Z1,ij , Z2,ij ) と書ける．ここで，• Tij = min(Xij , Dij , Cij )；最初に起こったイベントの時間，• δij = I(Tij = Xij )；増悪の指標（無増悪 = 0；増悪 = 1），• Tij∗ = min(Dij , Cij )；最後に起こったイベントの時間，∗• δij= I(Tij∗ = Dij )；死亡の指標（生存 = 0；死亡 = 1），• Z1,ij；増悪に関連する p1 次元共変量，• Z2,ij；死亡に関連する p2 次元共変量．
∗∗) は上記 A，B，C，
D の 4 つのケースを区別するための情報であり，(δij , δij) = (1, 1),
ペア (δij , δij∗（表 3）．例えば，ケース A は (δij , δij ) = (1, 1) に対応(1, 0), (0, 1), (0, 0) のいずれかを取りえる
し，この場合 Tij = Xij と Tij∗ = Dij が観測される．
上記のような構造をもつ生存期間データは「半競合リスクデータ」と呼ばれている（Fine et（terminal event）が増悪のような非終着イベン
al., 2001）．これは，死亡のような終着イベントの競合リスクになっていることに由来する．ここで「終着イベント」とト（non-terminal event）
は，それが起こった後の一切の情報が得られなくなるようなイベントを指す．死亡は明らかに
終着イベントである．増悪後には死亡時刻などの情報が得られるため，増悪は終着イベントでない．半競合リスクデータでは増悪だけでなく，より一般の中間イベント（非終着イベント）を考えることが出来る．半競合リスクデータの統計解析法の最新の動向は Peng（2019）を参照さ
れたい．

統計数理

158

第 68 巻 第 1 号 2020

表 3．患者がとりうる 4 つのケース（A〜D）．

5.2 Joint frailty-copula モデル
Rondeau et al.（2015）
は研究単位間の異質性を考慮した半競合リスクデータに基づくメタ分析の手法を提案している．彼女らが提案している Joint frailty モデルにおいて，Frailty 項
ui > 0 は Varη (ui ) = η のガンマ分布のモデル（4.1）に従うと仮定し，分散 η は研究単位間の異（ui ，Z1,ij ，Z2,ij を与えたとき質性の度合いを決めるパラメータである．Xij と Dij の条件付の）ハザード関数が次のような Joint frailty モデルに従うと仮定する．

rij (t|ui ) = ui r0 (t) exp(β1 Z1,ij )
λij (t|ui ) = uαi λ0 (t) exp(β2 Z2,ij ).

ここで β· は Z·,ij の効果を決める未知パラメータで，r0 (·) と λ0 (·) は未知の基準ハザード関数である．他の研究単位と比較して高リスクな研究単位は ui > 1，低リスクな研究単位は
0 < ui < 1 で与えられ，パラメータ α は増悪率と死亡率の異質性の差を表す．例えば，α = 0
は増悪率に異質性があるが死亡率には異質性がないモデル，α = 1 は増悪率と死亡率に同等の異質性があるモデルに対応する．Rondeau et al.（2015）のモデルにおいて重要な仮定は，(ui , Z1,ij , Z2,ij ) が与えられたもとで，Xij と Dij が条件付独立となることである．
Emura et al.（2017）は条件付独立の仮定を除いたコピュラのモデルである Joint frailty-copula
モデルを提案した．Cθ をコピュラとしたとき，Joint frailty-copula モデルは次の式で与えられる；⎧⎨ rij (t|ui ) = ui r0 (t) exp(β1 Z1,ij )

(5.1)

⎩


λij (t|ui ) = uαi λ0 (t) exp(β2 Z2,ij )Pr(Xij > x, Dij > y|ui ) = Cθ [SXij (x|ui ), SDij (y|ui )].

ここで θ は未知パラメータである．このモデルにおいて，Xij と Dij の生存関数は次のように書ける

SXij (x|ui ) = exp{−ui R0 (x) exp(β1 Z1,ij )}SDij (y|ui ) = exp{−uαi Λ0 (y) exp(β2 Z2,ij )}.
xy

ここで R0 (x) = 0 r0 (t)dt と Λ0 (y) = 0 λ0 (t)dt は累積基準ハザード関数である．
コピュラ Cθ は Xij と Dij の相関構造を定める．独立型 Cθ (v, w) = vw は Rondeau et al.
（2015）のモデルに帰着する．クレイトン型ではケンドール相関は τθ = θ/(θ + 2) であり，グン
ベル型ではケンドール相関は τθ = θ/(θ + 1) である．Burzykowski et al.（2001）のコピュラモデ
ルと同じ原理で，ケンドール相関は，無増悪期間が全生存期間の代替エンドポイントになりえるかどうかの指標の 1 つである．
∗
Joint frailty-copula モデル内の未知パラメータは観測されたデータ{(Tij , Tij∗ , δij , δij, Z1,ij ,Z2,ij ), i = 1, 2, . . . , G, j = 1, 2, . . . , Ni }を用いて最尤法で推定することが出来る．まず次のようコピュラを用いた生存時間解析

159

な記号を導入する；
Rij (t) = R0 (t) exp(β1 Z1,ij ),rij (t) = dRij (t)/dt = r0 (t) exp(β1 Z1,ij ),Λij (t) = Λ0 (t) exp(β2 Z2,ij ),λij (t) = dΛij (t)/dt = λ0 (t) exp(β2 Z2,ij ).

モデル（4.1）と（5.1）のもとで，対数尤度関数は次のように構成される；(α, η, θ, β1 , β2 , r0 , λ0 )
G

 N

i

=
i=1

∗
{δij log rij (Tij ) + δij
log λij (Tij∗ )}

j=1

 ∞
u

+ log

mi +αm∗
i

0

Ni


∗

ψθ [uRij (Tij ), uα Λij (Tij∗ )]δij ψθ∗ [uRij (Tij ), uα Λij (Tij∗ )]δij

j=1

× Θθ [uRij (Tij ), uα∗Λij (Tij∗ )]δij δij Dθ [uRij (Tij ), uα Λij (Tij∗ )]Ni
ここで mi = j=1
δij は i 番目の研究単位で増悪した患者数，m∗i =単位で死亡した患者数であり，Dθ [s, t] = Cθ [exp(−s), exp(−t)]，[1,0]ψθ [s, t] =[1,0]Dθ [s, t],Dθ [s, t]ψθ∗ [s, t] =

[0,1]

[0,1]

Dθ [s, t],Dθ [s, t]Θθ [s, t] =

Ni

j=1


fη (u)du ,∗δij
は i 番目の研究[1,1][s, t]Dθ [s, t]
,
[1,0]
[0,1]
Dθ [s, t]Dθ [s, t]

Dθ

[1,1]

Dθ [s, t] = −∂Dθ [s, t]/∂s, Dθ [s, t] = −∂Dθ [s, t]/∂t, Dθ= ∂ 2 Dθ [s, t]/∂s∂t はコピュラを微分して得られる. 対数尤度の導出は Emura et al.（2017, 2019b）を参照されたい.
例えばクレイトン型では，Dθ [s, t] = Aθ (s, t)−1/θ ，ψθ [s, t] = exp(θs)/Aθ (s, t)，ψθ∗ [s, t] =exp(θt)/Aθ (s, t)，Θθ [s, t] = 1 + θ とおき，ここで Aθ (s, t) = exp(θs) + exp(θt) − 1 である．これ
らを用いて対数尤度を計算する．
基準ハザード関数 r0 (·) と λ0 (·) のモデル化はいくつか提案されている．Emura et al.（2017）はスプライン関数でモデル化し，罰則付き最尤法で推定することを提案している．スプライン関数は単調増加（減少）型，凸型，凹型などのハザードを柔軟にモデル化できる（Shih andEmura, 2020）．Peng et al.（2018）と Peng（2019）は完全にノンパラメトリックなモデルを提はパラメトリックなモデルを考え，パラメータに事前分布を仮案している．Li et al.（2019）と Wu et al.（2020）はワイブル分布を用いている．
定したベイズ法を用いている．篠原（2019）
Schneider et al.（2020）
でもワイブル分布と区分的指数分布のモデルを考えているが，Rondeauと同様に，条件付独立性を仮定している．スプライン関数を用いる場合，推定値
et al.（2015）は罰則付き対数尤度のヘッセ行列から計算でき（Emura et al.,の標準誤差（SE）と信頼区間（CI），R パッケージ joint.Cox（Emura, 2019）で実行可能である．ノンパラメトリック
2017, 2019b）
最尤法を用いる場合，ヘッセ行列を利用すると SE を過小評価する傾向があるため，層化ブー．
トストラップ法が推奨されている（Peng, 2019）
5.3 卵巣癌データへの適用例
Ganzfried et al.（2013）の卵巣癌データを用いて，2 つのエンドポイントを同時に解析する方は IPD メタ分析を用いて，遺伝子 CXCL12 の発現量の法を説明する．Ganzfried et al.（2013）
全生存期間に与える影響が有意であることを示した．本稿では，Joint frailty-copula モデルの（無増悪期間と全生存期間）に与える影適用例として，CXCL12 遺伝子が 2 つのエンドポイント

160

統計数理

第 68 巻 第 1 号 2020

表 4．卵巣癌データ（Ganzfried et al., 2013; Emura et al., 2017）．研究単位は GEO（GeneExpression Omnibus）アクセス番号をあらわす．イベント発生率（%）は研究単位ごとに特定のイベント（増悪，死亡，打ち切り）の発生率を計算した．

表 5．卵巣癌データにおける無増悪期間と全生存期間を同時解析した結果．CXCL12 遺伝子
発現量の相対リスクを考えている．相対リスクが 1 以上とは，CXCL12 値が高いと増悪（死亡）リスクが上昇することを意味する．

響の解析事例を紹介する．
Ganzfried et al.（2013）のデータは多くの研究単位を含む IPD メタ分析用のデータであるが，多くの研究単位では 1 つのエンドポイントしか測定されていない．そのため，Emura et al.
（2017）で行われたように，無増悪期間と全生存期間の両方が測定されている 4 つの研究単位か．この取り出したデータは R パッケージ joint.Cox
ら 1,003 人の卵巣癌患者を取り出した（表 4）に収められている．Ganzfried et al.（2013）
に従い，CXCL12 遺伝子の発現量のみを共変量とした．
卵巣癌データを Joint frailty-copula のクレイトン型コピュラモデルに当てはめた結果を表 5に示す．CXCL12 遺伝子が全生存期間に与える影響は有意（P 値 <0.05）であり，その相対リス
は RR=1.18（95%CI: 1.08-1.29）である．また，CXCL12 遺伝子が無増悪期間
ク（Relative Risk）
であり，その影響は全生存期間のそれよりも大きい（RR=1.22,に与える影響も有意（P 値 <0.05）
95%CI: 1.13-1.32）．このことは，CXCL12 の発現量は卵巣癌患者の術後増悪を予測するバイ
オマーカーとなりうることを示している．これら相対リスクの値は，CXCL12 の発現量がSD 1 つ分増加したときのリスク増加を表す．コピュラのパラメータ推定値（θ̂ = 2.35, 95%CI:
は，増悪と死亡の間に中程度の相関があることを示す（τ̂ = 0.54, 95%CI: 0.38-0.70）．
1.90-2.90）
このことは，増悪で死亡をある程度予測できることを意味する．
独立型コピュラを適用した結果を表 5 に示す．この場合，増悪と死亡間の相関を考慮できな
いため，このモデルの妥当性は弱い．しかし，相対リスクの推定値は，クレイトン型のものと比べてほとんど変化がない．これは，相対リスクの推定値はコピュラの誤特定に関してある程コピュラを用いた生存時間解析

161

度の頑健性を持つことを示している．このような頑健性は Emura et al.（2017）のシミュレー
ション解析でも報告されている．
6.

動的予測と個別医療

6.1 生存期間の予測
生存時間解析における重要なテーマの 1 つが，生存期間の予測法の開発である．過去の患者の生存期間データをもとにして，新たな患者の生存期間を予測するためのモデル構築に関する手法は数多く提案されており，これは統計学の手法に加えて機械学習の手法も多い．例えば，乳癌（Gómez et al., 2016; Shukla et al., 2018），卵巣癌（Yoshihara et al., 2010; Emura et al.,
，膀胱癌（Guinney et al., 2017），その他の癌種（Kim et al., 2018; Choi et al., 2018）につ2018）
いて，患者の生存期間の予測式が提案されている．
正確な予測式を構築することは，患者本人やその主治医に対してより正確な人生設計や治療計画を考える礎となるため，個別医療（personalized medicine）の中核となる．個別医療の定義によると，それは遺伝子（分子）情報を基に健康増進を図る医療である（Schleidgen et al., 2013;．これは遺伝子（分子）情報の活用が患者の予後予測に重要な役割を果たすた
Hayes et al., 2014）
めである．
Waldron et al.（2014）
は卵巣癌患者の大規模な IPD メタ分析用データを用い，遺伝子を利用した 14 種類の予測法の予測精度を比較した．彼らは，12 種類の予測法が全生存期間の予測にに働くことを示したが，同時にその予測能力はそれほど高くない（C 統計量有意（P 値 < 0.05）= 0.56 − 0.60）
ことも指摘した．彼らはまた，これら生存期間予測法が臨床現場で実際に有効利用されるためには，さらなる改良や検証が必要であることも指摘している．
生存期間予測法を改良するためのアイデアは数多くあるが，以下では動的予測法とコピュラを IPD メタ分析用データに適用した予測手法を考える．
6.2 動的予測
では，通常の共変量（年齢，腫瘍径など）に加えて，中間イベン
動的予測（dynamic prediction）
ト（増悪など）の情報を予測因子に取り入れる（van Houwelingen and Putter, 2011）．治療開始
時点（例えば手術日や投与開始日）
を起点とする通常の予測の場合，手術日に測定される年齢や腫瘍径を共変量として扱い，手術日から時刻 w 後（例えば w = 5 年）に死亡する確率を予測す
る．動的予測では，手術後に経過観察した後にしか得られない増悪イベントの情報などを取りに設定し，時刻 t まで生存した患者が，時刻 t + w
入れるため，起点を術後（例えば術後 t 日目）
までに死亡する（または生存している）
確率を予測する（Thompson et al., 2007）．
いま D を全生存期間，X を無増悪期間，Z を遺伝子情報を含む共変量とする．共変量の値
Z は起点 t = 0（例えば手術日）に測定され，その後変化しないと仮定する（時間依存性共変量は考えない）．
に対して動的予測を考えている．患者 1，2 は時刻 t 以
図 5 はで 4 人の患者（患者 1〜患者 4）前に死亡しているので，予測の対象とはならない．患者 3，4 は時刻 t において生存している
ので，予測を行うことが出来る．患者 3 は時刻 t 以前に増悪しているので，増悪時刻 X = x で条件付けを行った死亡率を予測する．患者 4 は，時刻 t 以前には増悪していないので X > t で条件付けを行った死亡率を予測する．
動的予測を実践するにあたり，予測時刻 t は前もって決めておく必要がある．例えば，時刻t は術後 3 ヶ月ごとに取るなどと決めておく．この予測スケジュールは患者の状態によって変更してはならず，また患者の死亡や増悪に関連のある事象に依存してはならない．例えば，患

162

統計数理

第 68 巻 第 1 号 2020

図 5．増悪情報に基づく死亡率の動的予測．患者 1，2 は時刻 t 以前に死亡しているため，動的予測の対象とならない．患者 3，4 に対する動的予測は，時刻 t 以前に観測されている増悪情報に基づいて行われる．

者が何らかの不調を訴えて来院した時刻や，増悪した時刻を t に設定すると，妥当な予測は行えなくなる．
動的予測を定式化するために次のような 2 つの関数を利用する．
（Conditional Failure Function）（i）条件付死亡率関数
動的予測の基礎となるのが，次の 2 つの条件付死亡率関数
F (t, t + w|X > t, Z) = Pr(D ≤ t + w|D > t, X > t, Z),
F (t, t + w|X = x, Z) = Pr(D ≤ t + w|D > t, X = x, Z),x≤tである．2 つの式にある条件 {D ≥ t} は，予測の対象である患者が時刻 t において生存していることを前提としている．1 つ目の式にある条件 {X > t} は，時刻 t においてまだ増悪していないことを意味する．2 つ目の式にある条件 {X = x} は，時刻 t 以前に増悪したことを意味する．
この場合，増悪した時刻 X = x は既知なので予測式に組み入れることができる．条件付死亡率
関数は，これら条件下で，患者が時刻 t から t + w の間に死亡する確率である．通常の（動的でない）予測では予測時刻を t = 0 とし，w = 5 年などに設定して予測を行う．予測時刻を t = 0と設定すると，増悪の情報は不明なので，予測式は F (t, t + w|Z) = Pr(D ≤ t + w|D > t, Z) のみである．
増悪した患者は，そうでない患者に比べて，その死亡率が上昇する．したがって，直感的に不等式 F (t, t + w|X = x, Z) > F (t, t + w|X > t, Z) が全ての w > 0 に対して成立すると（数学的に成立考えられる（図 1）．X と D が正相関を持つ場合，この不等式は通常成立する（Z を与えたときの）条件付独立であれば，等式するかどうかはモデルに依存）．X と D がF (t, t + w|X = x, Z) = F (t, t + w|X > t, Z) が成立することは明らかである．

コピュラを用いた生存時間解析

163（ii）条件付ハザード関数（Conditional Hazard Function）
条件付ハザード関数は，条件付死亡率関数を用いて次のように表現できる．
λ(t|X = x, Z) = F (t, t + dw|X = x, Z)/dw,
λ(t|X > x, Z) = F (t, t + dw|X > x, Z)/dw.

これらの式により，条件付ハザード関数は，上記の条件付死亡率関数において w が微小な値dw をとる場合と解釈できる．2 つの条件付ハザード関数の比をとると，X と D の関連を表すとなる．
指標である次のようなオッズ比関数（3.3 節）Pr(X = x, D = t|Z) Pr(X > x, D > t|Z)λ(t|X = x, Z)=,
λ(t|X > x, Z)
Pr(X = x, D > t|Z) Pr(X > x, D = t|Z)

t > 0,

x > 0.

オッズ比が 1 より大きいときは正相関を表し，オッズ比が 1 より小さいときは負相関を表す．
（Z を与えたときの）条件付独立である．
オッズ比が 1 のときは X と D はClayton（1978）
が提案した次のセミパラメトリックモデル
1

Pr(X > x, D > t|Z) = [SX (x|Z)−θ + SD (t|Z)−θ − 1]− θ ,t > 0,

x > 0,

θ>0を考えると．このモデルの下では，2 つの条件付ハザード関数の関係をオッズ比定数 (θ + 1) を用いてλ(t|X = x, Z) = (θ + 1)λ(t|X > x, Z), t > 0, x > 0.と表すことができる．これは，患者の増悪イベント {X = x} があったときと，そうでないとき {X > x} の相対リスクが (θ + 1) となるクレイトン型（Clayton, 1978）コピュラモデルである．Day et al.（1997）
はクレイトン型モデルにランドマーク解析を使った動的予測（3.3 節を参照）
を提案している．
動的予測において，条件付死亡率関数，条件付ハザード関数のどちらも有用な予測指標である．時刻 t を予測時刻と呼び，これは臨床医が患者に対して予測を始める時刻である．これら予測指標をデータから推定する場合，通常何らかのモデルを利用する．van Houwelingen andはランドマーク解析法を採用している．これは予測時刻 t 以前に収集されたデー
Putter（2011）
タを基に条件付 Cox 回帰法を当てはめて予測モデルを構築するものである．ランドマーク解
析法では，t を複数選び，複数の条件付 Cox モデルを構築したときにモデル間の一致性が成立しないなどのいくつかの諸問題がある．詳細については van Houwelingen and Putter（2011）とSuresh（2018）に詳しく議論されている．
ランドマーク解析法の諸問題の解決策として，結合モデル（Joint model）を用いた動的予測式が提案されている（Proust-Lima et al., 2014; Rondeau et al., 2017; Sène et al., 2016 など）．これらのモデルにおいては，Frailty 項を用いて応答変数間の相関を考慮している．Frailty 項の替
わりに，コピュラで相関を考慮するモデルは Emura et al.（2018），Suresh et al.（2019）などを参照されたい．
6.3 Joint frailty-copula モデルを利用した動的予測
以下では，Frailty 項とコピュラで構成した Joint frailty-copula モデルによる動的予測法を紹介する（Emura et al., 2018）．この手法では，全生存期間と無増悪期間の相関を解析する際に用いる IPD メタ分析データを利用する．このため，Frailty 項は相関を考慮するというよりも，研究単位間の異質性をモデル化するために導入されている．
無増悪期間の生存関数を SX (x|u) = Pr(X > x|u, Z)，全生存期間の生存関数を SD (y|u) =

統計数理

164

第 68 巻 第 1 号 2020

Pr(D > y|u, Z) で定義し，ここで u は共変量 Z で説明しきれない Frailty 項である．u の分布
は，式（4.1）で与えられるように，平均 1，分散 η のガンマ分布であるとする．Joint frailty-copula
モデルの仮定をおくと，Pr(X > x, D > y|u, Z ) = Cθ [SX (x|u), SD (y|u)],となり，Cθ (v, w) はコピュラ，θ は X と D の相関を表すパラメータである．このモデルの下で，条件付死亡率関数は次のように表せる．
（I）時刻 t 以前に患者が増悪していない場合
X > t で条件付けを行い，条件付死亡率関数はF (t, t + w|X > t, Z) = Pr(D ≤ t + w|D > t, X > t, Z)∞(Cθ [SX (t|u), SD (t|u)] − Cθ [SX (t|u), SD (t + w|u)])fη (u)du= 0∞
Cθ [SX (t|u), SD (t|u)]fη (u)du0となる．ただし fη (u) は式（4.1）
で与えられるガンマ分布の確率密度関数である．また，条件付
ハザード関数は[0,1]

λ(t|X > t, Z, u) = λD (t|u)

SD (t|u)Cθ [SX (t|u), SD (t|u)]Cθ [SX (t|u), SD (t|u)]となる．ここで λD (t|u) = −∂ log SD (t|u)/∂t，Cθ[0,1] (v, w) = ∂Cθ (v, w)/∂w である．
（II）時刻 t 以前に患者が増悪している場合
増悪時刻 X = x で条件付けを行い，条件付死亡率関数はF (t, t + w|X = x, Z) = Pr(D ≤ t + w|D > t, X = x, Z)
[1,0]
[1,0]
∞
(Cθ [SX (x|u), SD (t|u)] − Cθ [SX (x|u), SD (t + w|u)])λX (x|u)SX (x|u)fη (u)du
= 0
[1,0]
∞
Cθ [SX (x|u), SD (t|u)]λX (x|u)SX (x|u)fη (u)du0となる．ここで λX (x|u) = −∂ log SX (x|u)/∂x である．また，条件付ハザード関数はλ(t|X = x, Z, u) = λD (t|u)[1,1][SX (x|u), SD (t|u)][1,0]Cθ [SX (x|u), SD (t|u)]

SD (t|u)Cθ

となる．ここで Cθ[1,0] (v, w) = ∂Cθ (v, w)/∂v ，Cθ[1,1] (v, w) = ∂ 2 Cθ (v, w)/∂v∂w である．
を参照されたい．IPD メタ分析用データなどを用上式の導出は Emura et al.（2018, 2019b）
いて条件付死亡率関数や条件付ハザード関数内のパラメータは推定することができる（詳しく．
は Emura et al., 2018, 2019b を参照）
条件付死亡率関数は Frailty 項 u で積分している．これは予測すべき患者の Frailty 項 u が観測不能なため，u の分布で平均化していることを意味する．積分は数値的に計算する．簡便な方法であるが，積分する代わりに平均値 u = 1 で置き換える方法もある．
条件付ハザード関数は Frailty 項 u で条件付けを行っている．u は未知なので積分する方が望ましいが，条件付けを行うことにより，次のような数学的性質が利用できる．
λ(t|X = t, Z, u)= Rθ [SX (t|u), SD (t|u)].
λ(t|X > t, Z, u)

ここで Rθ は 3.3 節で既出のオッズ比関数である．クレイトン型コピュラでは Rθ (v, w) = 1 + θが {X = t} の {X > t} に対する相対リスクとなる．先に λ(t|X > t, Z, u) を計算すれば，それに 1 + θ をかけることにより λ(t|X = t, Z, u) が得られる．予測に使用する Frailty 項は u = 1

コピュラを用いた生存時間解析

165

とする．
動的予測を臨床現場で実践する場合，まだ死亡していない新たな患者を対象にする．ここで「新たな」とは，すでに死亡・生存の状態が確定しているデータ内の患者と区別するための語である．現実にそのような患者を有していなくても，仮想的な患者を想定し，どのように予測のプロセスが実行されるのかを考察することは重要である．
まず患者が有している共変量情報を Z = (Z1 , Z2 , CC1 , CC2 ) とし，これらは：• Z1；増悪に関連する p1 次元の臨床的共変量，• Z2；全生存に関連する p2 次元の臨床的共変量，• CC1 = w1 V1 + · · · + wq1 Vq1；増悪に関連する遺伝子発現量の複合共変量，• CC2 = 1 W1 + · · · + q2 Wq2；全生存に関連する遺伝子発現量の複合共変量，と定義する．ここで (V1 , . . . , Vq1 ) は q1 個の増悪に関連する遺伝子発現量，(W1 , . . . , Wq2 ) は q2個の全生存に関連する遺伝子発現量である．遺伝子の個数 q（と重み w（1 または q2 ）k またはk）はデータから計算される（Emura et al., 2018, 2019b）．全ての遺伝子発現量は，遺伝子ごとに（compound covariate）の基本
平均 0，標準偏差 1 で標準化されていると仮定する．複合共変量，Emura et al.（2012, 2019a, 2019b）を参照されたい．
的な考え方については Matsui（2006）
データを Joint frailty-copula モデルに当てはめることにより，無増悪期間と全生存期間の生存関数が次のように推定される．







CC1 − μ̂1
,
σ̂1



CC2 − μ̂2
ŜD (t|u) = exp −uα Λ̂0 (t) exp β̂2 Z2 + γ̂2.
σ̂2ŜX (t|u) = exp −uR̂0 (t) exp β̂1 Z1 + γ̂1tt

こ こ で R̂0 (t) = 0 r̂0 (x)dx と Λ̂0 (t) = 0 λ̂0 (x)dx は 推 定 さ れ た 累 積 基 準 ハ ザ ー ド 関 数 ，(θ̂, η̂, β̂1 , β̂2 , γ̂1 , γ̂2 , r̂0 , λ̂0 ) は推定値されたパラメータの値である．μ̂· は CC· の平均，σ̂· はCC· の標準偏差である．
これら推定値を条件付死亡率関数や条件付ハザード関数の未知パラメータに代入することにより，新たな患者への動的予測が可能となる．例えば，条件付死亡率関数は次のように計算される．
∞

F̂ (t, t + w|X > t, Z) =0

(Cθ̂ [ŜX (t|u), ŜD (t|u)] − Cθ̂ [ŜX (t|u), ŜD (t + w|u)])fη̂ (u)du∞0Cθ̂ [ŜX (t|u), ŜD (t|u)]fη̂ (u)du,F̂ (t, t + w|X = x, Z)
∞

=

0

[1,0]

(Cθ̂

[1,0]

[ŜX (x|u), ŜD (t|u)] − Cθ∞0[ŜX (x|u), SD (t + w|u)])λ̂X (x|u)ŜX (x|u)fη̂ (u)du[1,0]Cθ̂ [ŜX (x|u), ŜD (t|u)]λ̂X (x|u)ŜX (x|u)fη̂ (u)du.

6.4 卵巣癌患者への適用例
Emura et al.（2018）
は Ganzfried et al.（2013）の 912 人の卵巣癌患者をもとに，次のような複合共変量を得た．
• 増悪に関連する 158 個の遺伝子（P 値 < 0.001）；
CC1 = (0.249 × CXCL12 ) + (0.235 × TIMP2 ) + (0.222 × PDPN ) + · · · + (−0.152 × MMP12 ).
• 死亡に関連する 128 個の遺伝子（P 値 < 0.001）；

166

統計数理

第 68 巻 第 1 号 2020

CC2 = (0.237×NCOA3 )+(0.223×TEAD1 )+(0.263×YWHAB)+· · ·+(−0.157×KCNH4 ).

ここで遺伝子発現量は遺伝子名（CXCL12 など）で表し，全患者で平均 0，標準偏差 1 となるよう標準化してある．CC1 の平均は μ̂1 = 0.338，標準偏差は σ̂1 = 10.468 であり，CC2 の平均はμ̂2 = 0.222，標準偏差は σ̂2 = 7.894 であった．P 値が 0.001 より小さい遺伝子を選択する基準
はマイクロアレイ解析でしばしば利用される（Matsui et al., 2012）．
データを Joint frailty-copula モデルに当てはめた結果，Pr(X > x, D > y|u ) = Cθ̂ [ŜX (x|u), ŜD (y|u)] = [ŜX (x|u)−θ̂ + ŜD (y|u)−θ̂ − 1]−1/θ̂ ,となり，ここで θ̂ = 1.9（ケンドール相関 τ̂ = 0.49），



ŜX (x|u) = exp −uR̂0 (t) exp γ̂1CC1 − μ̂1
σ̂1


,



CC2 − μ̂2,σ̂2となった．ここで Z2 は（0 または 1）をとる術後腫瘍径（≤ 1 cm or > 1 cm）である．腫瘍径が無増悪期間に与える影響は有意でなかった（P 値>0.05）ので ŜX (x|u) のモデルから除いてある．
回帰推定値は γ̂1 = 0.39，β̂2 = 0.16，γ̂2 = 0.44 であった．基準ハザード関数の推定値は，M.(t)を M スプライン基底関数（Emura et al., 2017, 2019b）とすると，ŜD (y|u) = exp −Λ̂0 (t) exp β̂2 Z2 + γ̂2r̂0 (t) = dR̂0 (t)/dt = 0.85 × M1 (t) + 2.14 × M2 (t) + 0 × M3 (t) + 0.07 × M4 (t) + 0 × M5 (t),λ̂0 (t) = dΛ̂0 (t)/dt = 0.17 × M1 (t) + 1.05 × M2 (t) + 1.24 × M3 (t) + 0.27 × M4 (t) + 0 × M5 (t),t ∈ [0, 6420] であり，6420（日）
は全患者の最長追跡日数である．異質性分散の推定値は Var(u)= η̂ = 0.04 であった．ステージなどの他の共変量もいくつか測定されているが，これらを追加してもモデルの当てはまりは改善しなかった．
次に，以下のような仮想的な二人の卵巣癌患者を考え，どのように予測のプロセスが実行されるのかを説明する．
• 患者 1；高リスク遺伝子（CC1 = 10，CC2 = 10），腫瘍径 > 1 cm（Z2 = 1）．
• 患者 2；低リスク遺伝子（CC1 = −10，CC2 = −10），腫瘍径 ≤ 1 cm（Z2 = 0）．

例えば患者 2 に対して条件付死亡率関数 F̂ (t, t + w|X = x, Z) を t = 1000，t + w < 6420，
x = 600 の設定で計算するときは，R パッケージ joint.Cox を利用し，「付録」にあるコードを実行すれば良い．CC を計算するには，患者の 158 個の遺伝子を調べ，重み付け和を計算する必要がある．ここではそのプロセスを省略した．
図 6 は二人の患者の条件付死亡率の推移を表す．t = 500 日目を起点に予測を始めると，患（図 6 の左パネル）．これは患者 1 が患者 2 より高リス
者 1 は患者 2 よりも高い死亡率を示す
ク遺伝子をもち，かつ腫瘍径が大きいことに由来する．図 6 の左パネルにおいては予測時刻t = 500 日目において両患者ともに増悪がないと仮定した．時刻を進め，t = 1000 日目を起点（x = 600 日目で）と仮定すると，患者 2 の死亡率が急激に上とし，患者 2 にのみ増悪があった（図 6 の右パネル）．このことは，増悪の死亡に与える影響が，昇し，患者 1 の死亡率を上回る
遺伝子と腫瘍径の死亡に与える影響よりも強いことを示している．
図 7 は二人の卵巣癌患者の条件付ハザード関数 λ̂(t|X = t, Z, u) と λ̂(t|X > t, Z, u) を u = 1で計算したものを示している．患者 1 は危険因子（高リスク遺伝子と腫瘍径>1cm）を持つため，患者 2 よりも高いハザード値を示していることがわかる．他方，患者 2 は時間とともになだらかに減少するハザード値を示している．同一の患者で，増悪した場合としなかった場合を比較

コピュラを用いた生存時間解析

167

図 6．二人の卵巣癌患者の条件付死亡率の推移．

図 7．二人の卵巣癌患者の条件付ハザードの推移．

するときに，これら 2 つのハザード間に比例関係が存在し，λ̂(t|X = t, Z, u) = (1 + θ̂)λ̂(t|X > t, Z, u) = 2.9λ̂(t|X > t, Z, u)と表せることに注意する．これはクレイトン型コピュラの仮定から導かれるものである．したがって，増悪した場合は，増悪しなかった場合と比較してハザードが約 3 倍になることがわかる．逆に増悪しなければ，患者 2 の死亡リスクはかなり低下することがわかる．また，増悪しなかった場合でも，患者 1 のように危険因子（高リスク遺伝子と腫瘍径 > 1 cm）を持つ患者は，将来的に死亡リスクが上昇する可能性があることがわかる．
7.

今後の課題

コピュラの型は無数にあるため，研究者はその目的に応じて最も適したコピュラを選択することができる．しかし実際には，生存期間において変数間の相関をモデリングするためにはクレイトン型のコピュラが最も広く利用されてきた．本稿で示したようにその数学的に簡潔な関数形，ケンドール相関やオッズ比の特徴からその利便性・有用性は理解されるであろう．また

統計数理

168

第 68 巻 第 1 号 2020

相関パラメータが相対リスクとして解釈できる医学上の重要性も強調した．このような理由から，本稿でも主にクレイトン型コピュラを用いて解説してきたが，これは母集団がクレイトン
型コピュラに従うことを必ずしも意味していない．母集団の相関構造をよりよく近似するコピュラを用いたい場合は，クレイトン型，グンベル型，フランク型など特徴が異なるいくつかのコピュラ間でモデル選択を行うことが望ましい．しかしながら，尤度や予測モデルの計算が複雑になることから，現時点では 5 節，6 節の手法は R パッケージ joint.Cox においてクレイトン型のみ利用可能となっている．
6 節においては，患者の生存期間の予測にコピュラをどのように利用するのかを，動的予測法の枠組みで解説した．患者の生存期間を正確に予測することは，個別医療の観点からも重要な課題であるため，現場の臨床医や医師が利用できる状態にすることが重要である．1 つの取り組みとして，臨床医向けに動的予測の利用法ウェブ・アプリケーションを作成する手法を解説した論文を執筆している（Emura et al., 2020）．
国内における上記動的予測の手法は，間質性肺炎合併肺癌患者へ応用された例がある（高見．またこの動的予測の手法をワイブル分布で行うことにより，長期予測法を実現する
澤, 2018）．
方法も提案されている（篠原, 2019）謝辞

本特集号の企画者である吉羽要直先生と加藤昇吾先生には，執筆の機会を頂いたことに心より感謝申し上げます．また改稿に当たり，有益なコメントをいただいた査読者と編集委員の方に感謝申し上げます．最後に，図 1 をご提供頂いた篠原紗也華さんに感謝を申し上げます．

参

考

文

献

Alonso, A. and Molenberghs, G. (2008). Evaluating time to cancer recurrence as a surrogate marker
for survival from an information theory perspective, Statistical Methods in Medical Research,
17(5), 497-504.
Barthel, N., Geerdens, C., Killiches, M., Janssen, P. and Czado, C. (2018). Vine copula based likelihood
estimation of dependence patterns in multivariate event time data, Computational Statistics &Data Analysis, 117, 109-127.
Burzykowski, T., Molenberghs, G., Buyse, M., Geys, H. and Renard, D. (2001). Validation of surrogate
end points in multiple randomized clinical trials with failure time end points, Applied Statistics,
50(4), 405-422.
Burzykowski, T., Molenberghs, G. and Buyse, M. (eds.) (2005). The Evaluation of Surrogate Endpoints, Springer, New York.
Buyse, M., Burzykowski, T., Michiels, S. and Carroll, K. (2008). Individual-and trial-level surrogacy
in colorectal cancer, Statistical Methods in Medical Research, 17, 467-475.
Chaieb, L., Rivest, L. P. and Abdous, B. (2006). Estimating survival under a dependent truncation,
Biometrika, 93(3), 655-669.
Chen, Y. H. (2012). Maximum likelihood analysis of semicompeting risks data with semiparametric
regression models, Lifetime Data Analysis, 18, 36-57.
Choi, J., Oh, I., Seo, S. and Ahn, J. (2018). G2Vec: Distributed gene representations for identiﬁcationof cancer prognostic genes, Scientific Reports, 8(1), 13729.
Clayton, D. G. (1978). A model for association in bivariate life tables and its application in epidemiological studies of familial tendency in chronic disease incidence, Biometrika, 65(1), 141-151.

コピュラを用いた生存時間解析

169

Cox, D. R. (1972). Regression models and life-tables (with discussion), Journal of the Royal StatisticalSociety, Series B, 34, 187-220.
Day, R., Bryant, J. and Lefkopoulou, M. (1997). Adaptation of bivariate frailty models for prediction,with application to biological markers as prognostic indicators, Biometrika, 84(1), 45-56.
de Uña-Álvarez, J. and Veraverbeke, N. (2013). Generalized copula-graphic estimator, TEST, 22(2),343-360.
de Uña-Álvarez, J. and Veraverbeke, N. (2017). Copula-graphic estimation with left-truncated andright-censored data, Statistics, 51(2), 387-403.
Duchateau, L., Janssen, P., Lindsey, P., Legrand, C., Nguti, R., et al. (2002). The shared frailty modeland the power for heterogeneity tests in multicenter trials, Computational Statistics & Data
Analysis, 40(3), 603-620.
Eisenhauer, E. A., Therasse, P., Bogaerts, J., Schwartz, L. H., Sargent, D., et al. (2009). New response
evaluation criteria in solid tumours: Revised RECIST guideline (version 1.1), European Journalof Cancer, 45(2), 228-247.
Emura, T. (2019). joint.Cox: Joint frailty-copula models for tumour progression and death in metaanalysis, CRAN.
Emura, T. and Chen, Y. H. (2016). Gene selection for survival data under dependent censoring, a
copula-based approach, Statistical Methods in Medical Research, 25(6), 2840-2857.
Emura, T. and Chen, Y. H. (2018). Analysis of Survival Data with Dependent Censoring, Copula-based
Approaches, JSS Research Series in Statistics, Springer, Singapore.
Emura, T. and Michimae, H. (2017). A copula-based inference to piecewise exponential models under
dependent censoring, with application to time to metamorphosis of salamander larvae, Environmental and Ecological Statistics, 24(1), 151-173.
Emura, T. and Murotani, K. (2015). An algorithm for estimating survival under a copula-based dependent truncation model, TEST, 24(4), 734-751.
Emura, T. and Pan, C. H. (2020). Parametric likelihood inference and goodness-of-ﬁt for dependently
left-truncated data, a copula-based approach, Statistical Papers, 61(1), 479-501.
Emura, T. and Wang, W. (2010). Testing quasi-independence for truncation data, Journal of Multivariate Analysis, 101, 223-239.
Emura, T. and Wang, W. (2012). Nonparametric maximum likelihood estimation for dependent truncation data based on copulas, Journal of Multivariate Analysis, 110, 171-188.
Emura, T., Lin, C. W. and Wang, W. (2010). A goodness-of-ﬁt test for Archimedean copula models in
the presence of right censoring, Computational Statistics & Data Analysis, 54, 3033-3043.
Emura, T., Wang, W. and Hung, H. N. (2011). Semi-parametric inference for copula models for dependently truncated data, Statistica Sinica, 21, 349-367.
Emura T., Chen, Y. H. and Chen, H. Y. (2012). Survival prediction based on compound covariate underCox proportional hazard models, PLoS One, 7(10), e47627, DOI:10.1371/journal.pone.0047627.
Emura, T., Nakatochi, M., Murotani, K. and Rondeau, V. (2017). A joint frailty-copula model between tumour progression and death for meta-analysis, Statistical Methods in Medical Research,
26(6), 2649-2666.
Emura, T., Nakatochi, M., Matsui, S. and Michimae, H. and Rondeau, V. (2018). Personalized dynamic prediction of death according to tumour progression and high-dimensional genetic factors:
Meta-analysis with a joint model, Statistical Methods in Medical Research, 27(9), 2842-2858.
Emura, T., Matsui, S. and Chen, H. Y. (2019a). compound.Cox: Univariate feature selection andcompound covariate for predicting survival, Computer Methods and Programs in Biomedicine,
168, 21-37.
Emura, T., Matsui, S. and Rondeau, V. (2019b). Survival Analysis with Correlated Endpoints, Joint

170

統計数理

第 68 巻 第 1 号 2020

Frailty-Copula Models, JSS Research Series in Statistics, Springer, Singapore.
Emura, T., Shih, J. H., Ha, I. D. and Wilke, R. A. (2019c). Comparison of the marginal hazard modeland the sub-distribution hazard model for competing risks under an assumed copula, Statistical
Methods in Medical Research, DOI: 10.1177/0962280219892295.
Emura, T., Michimae, H. and Matsui, S. (2020). Building web applications for personalized risk prediction using R packages joint.Cox and Shiny (submitted).
Fine, J. P., Jiang, H. and Chappell, R. (2001). On semi-competing risks data, Biometrika, 88, 907-920.
Ganzfried, B. F., Riester, M., Haibe-Kains, B., Risch, T., Tyekucheva, S., et al. (2013). Curated
ovarian data: Clinically annotated data for the ovarian cancer transcriptome, Database, Article
ID bat013, DOI:10.1093/database/bat013.
Gómez, I., Ribelles, N., Franco, L., Alba, E. and Jerez, J. M. (2016). Supervised discretization can discover risk groups in cancer survival analysis, Computer Methods and Programs in Biomedicine,
136, 11-19.
Guinney, J., Wang, T., Laajala, T. D., Winner, K. K., Bare, J., et al. (2017). Prediction of overall
survival for patients with metastatic castration-resistant prostate cancer: Development of aprognostic model through a crowdsourced challenge with open clinical trial data, The Lancet
Oncology, 18(1), 132-142.
Gumbel, E. J. (1960). Distributions de valeurs extremes en plusieurs dimensions, Publications del’Institut de statistique de l’Université de Paris, 9, 171-173.
Hayes, D. F., Markus, H. S., Leslie, R. D. and Topol, E. J. (2014). Personalized medicine: Risk
prediction, targeted therapies and mobile health technology, BMC Medicine, 12(1), 37.
Hirsch, K. and Wienke, A. (2012). Software for semiparametric shared gamma and log-normal frailty
models: An overview, Computer Methods and Programs in Biomedicine, 107(3), 582-597.
Kim, M., Oh, I. and Ahn, J. (2018). An improved method for prediction of cancer prognosis by network
learning, Genes, 9, 478.
Klein, J. P., Van Houwelingen, H. C., Ibrahim, J. G. and Scheike, T. H. (eds.) (2014). Handbook ofSurvival Analysis, CRC Press, New York.
Le Tourneau, C., Michiels, S., Gan, H. K. and Siu, L. L. (2009). Reporting of time-to-event end points
and tracking of failures in randomized trials of radiotherapy with or without any concomitant anticancer agent for locally advanced head and neck cancer, Journal of Clinical Oncology,27(35), 5965-5971.
Li, Z., Chinchilli, V. M. and Wang, M. (2019). A Bayesian joint model of recurrent events and aterminal event, Biometrical Journal, 61(1), 187-202.
Li, Z., Chinchilli, V. M. and Wang, M. (2020). A time-varying Bayesian joint hierarchical copula
model for analysing recurrent events and a terminal event: An application to the Cardiovascular Health Study, Journal of the Royal Statistical Society: Series C (Applied Statistics), 69(1),151-166.
Ling, M., Hu, T. and Sun, J. (2016). Cox regression analysis of dependent interval-censored failure
time data, Computational Statistics & Data Analysis, 103, 79-90.
Lo, S. M. and Wilke, R. A. (2010). A copula model for dependent competing risks, Journal of theRoyal Statistical Society, Series C, 59(2), 359-376.
Matsui, S. (2006). Predicting survival outcomes using subsets of signiﬁcant genes in prognostic marker
studies with microarrays, BMC Bioinformatics, 7, 156.
Matsui, S., Simon, R., Qu, P., Shaughnessy, J. D., Barlogie, B. and Crowley, J. (2012). Developing andvalidating continuous genomic signatures in randomized clinical trials for predictive medicine,
Clinical Cancer Research, 18(21), 6065-6073.
Michiels, S., Le Maître, A., Buyse, M., Burzykowski, T., Maillard, E., et al. (2009). Surrogate end-

コピュラを用いた生存時間解析

171

points for overall survival in locally advanced head and neck cancer: Meta-analyses of individual
patient data, Lancet Oncology, 10(4), 341-350.
Michiels, S., Pugliano, L., Marguet, S., Grun, D., Barinoﬀ, J., et al. (2016). Progression-free survival as surrogate end point for overall survival in clinical trials of HER2-targeted agents in
HER2-positive metastatic breast cancer, Annals of Oncology, 27(6), 1029-1034.
Moradian, H., Larocque, D. and Bellavance, F. (2019). Survival forests for data with dependent censoring, Statistical Methods in Medical Research, 28(2), 445-461.
Nelsen, R. B. (2006). An Introduction to Copulas, 2nd ed., Springer, New York.
Oakes, D. (1989). Bivariate survival models induced by frailties, Journal of the American Statistical
Association, 84, 487-493.
Oba, K., Paoletti, X., Alberts, S., Bang, Y. J., Benedetti, J., et al. (2013). Disease-free survival as asurrogate for overall survival in adjuvant trials of gastric cancer: A meta-analysis, Journal of
the National Cancer Institute, 105(21), 1600-1607.
Pazdur, R. (2008). Endpoints for assessing drug activity in clinical trials, The Oncologist, 13, 19-21.
Peng, M. (2019). Analysis of Complex Survival Data Subject to Semi-competing Risks, Doctoral Thesis,
Nanyang Technological University, Singapore.
Peng, M. and Xiang, L. (2019). Joint regression analysis for survival data in the presence of two setsof semi-competing risks, Biometrical Journal, DOI: 10.1002/bimj.201800137.
Peng, M., Xiang, L. and Wang, S. (2018). Semiparametric regression analysis of clustered survival datawith semicompeting risks, Computational Statistics & Data Analysis, 124, 53-70.
Piedbois, P. and Croswell, M. J. (2008). Surrogate endpoints for overall survival in advanced colorectal
cancer: A clinician’s perspective, Statistical Methods in Medical Research, 17(5), 519-527.
Proust-Lima, C., Séne, M., Taylor, J. M. and Jacqmin-Gadda, H. (2014). Joint latent class models forlongitudinal and time-to-event data: A review, Statistical Methods in Medical Research, 23(1),74-90.
Riley, R. D. (2009). Multivariate meta-analysis: The eﬀect of ignoring within-study correlation, Journalof the Royal Statistical Society, Series A, 172, 789-811.
Rondeau, V. and Gonzalez, J. R. (2005). frailtypack: A computer program for the analysis of correlated
failure time data using penalized likelihood estimation, Computer Methods and Programs in
Biomedicine, 80(2), 154-164.
Rondeau, V., Pignon, J. P. and Michiels, S. (2015). A joint model for dependence between clustered
times to tumour progression and deaths: A meta-analysis of chemotherapy in head and neckcancer, Statistical Methods in Medical Research, 24(6), 711-729.
Rondeau, V., Mauguen, A., Laurent, A., Berr, C. and Helmer, C. (2017). Dynamic prediction models
for clustered and interval-censored outcomes: Investigating the intra-couple correlation in therisk of dementia, Statistical Methods in Medical Research, 26(5), 2168-2183.
Rotolo, F., Legrand, C. and Van Keilegom, I. (2013). A simulation procedure based on copulas togenerate clustered multi-state survival data, Computer Methods and Programs in Biomedicine,
109(3), 305-312.
Rotolo, F., Paoletti, X. and Michiels, S. (2018). surrosurv: An R package for the evaluation of failure
time surrogate endpoints in individual patient data meta-analyses of randomized clinical trials,
Computer Methods and Programs in Biomedicine, 155, 189-198.
Schleidgen, S., Klingler, C., Bertram, T., Rogowski, W. H. and Marckmann, G. (2013). What is personalized medicine: Sharpening a vague term based on a systematic literature review, BMC
Medical Ethics, 14(1), 55.
Schneider, S., Demarqui, F. N., Colosimo, E. A. and Mayrink, V. D. (2020). An approach to modelclustered survival data with dependent censoring, Biometrical Journal, 62(1), 157-174.

172

統計数理

第 68 巻 第 1 号 2020

Sène, M., Taylor, J. M., Dignam, J. J., Jacqmin-Gadda, H. and Proust-Lima, C. (2016). Individualized dynamic prediction of prostate cancer recurrence with and without the initiation of asecond treatment: Development and validation, Statistical Methods in Medical Research, 25(6),2972-2991.
Sherrill, B., Amonkar, M., Wu, Y., Hirst, C., Stein, S, et al. (2008). Relationship between eﬀects ontime-to-disease progression and overall survival in studies of metastatic breast cancer, British
Journal of Cancer, 99, 1542-1548.
Shi, Q. and Sargent, D. J. (2009). Meta-analysis for the evaluation of surrogate endpoints in cancer
clinical trials, International Journal of Clinical Oncology, 14(2), 102-111.
Shih, J. H. and Emura, T. (2018). Likelihood-based inference for bivariate latent failure time models with competing risks under the generalized FGM copula, Computational Statistics, 33(3),1293-1323.
Shih, J. H. and Emura, T. (2020). Penalized Cox regression with a ﬁve-parameter spline model (submitted).
Shih, J. H., Lee, W., Sun, L. H. and Emura, T. (2019). Fitting competing risks data to bivariate Paretomodels, Communications in Statistics — Theory and Methods, 48(5), 1193-1220.
篠原紗也華 (2019). Dynamic prediction using joint frailty-copula model: Long-term prediction by
Weibull distribution, 薬学学士論文, 北里大学薬学部薬学科臨床統計学教室.
Shukla, N., Hagenbuchner, M., Win, K. T. and Yang, J. (2018). Breast cancer data analysis forsurvivability studies and prediction, Computer Methods and Programs in Biomedicine, 155,
199-208.
Sklar, A. (1959). Fonctions de répartition à n dimensions et leurs marges, Publications de l’Institut deStatistique de L’Université de Paris, 8, 229-231.
Sofeu, C., Emura, T. and Rondeau, V. (2019). One-step validation method for surrogate endpoints in
multiple randomized cancer clinical trials with failure-time endpoints, Statistics in Medicine,
38(16), 2928-2942.
Sofeu, C., Emura, T. and Rondeau, V. (2020). A joint frailty-copula model for meta-analytic validationof failure time surrogate endpoints in clinical trials (in revision).
Suresh, K. (2018). Dynamic Methods for the Prediction of Survival Outcomes using Longitudinal
Biomarkers, Ph.D Dissertation, University of Michigan Library.
Suresh, K., Taylor, J. M., Tsodikov, A. (2019). A Gaussian copula approach for dynamic prediction ofsurvival with a longitudinal biomarker, Biostatistics, DOI: 10.1093/biostatistics/kxz049.
高見澤勇 (2018). 化学療法が対象となる間質性肺炎合併肺がん患者の死亡率の動的予測：Joint frailtycopula モデルの臨床データへの適用, RIMS Kôkyûroku, No. 2091, 31-37.
Thompson, R. H., Leibovich, B. C., Lohse, C. M., Cheville, J. C., Zincke, H., et al. (2007). Dynamic outcome prediction in patients with clear cell renal cell carcinoma treated with radical
nephrectomy: The D-SSIGN score, The Journal of Urology, 177(2), 477-480.
Yoshihara, K., Tajima, A., Yahata, T., Kodama, S., Fujiwara, H., et al. (2010). Gene expression
proﬁle for predicting survival in advanced-stage serous ovarian cancer across two independent
datasets, PloS One, 5(3), e9615.
van Houwelingen, H. C. and Putter, H. (2011). Dynamic Prediction in Clinical Survival Analysis, CRC
Press, New York.
Waldron, L., Haibe-Kains, B., Culhane, A. C., Riester, M., Ding, J., et al. (2014). Comparative metaanalysis of prognostic gene signatures for late-stage ovarian cancer, Journal of the National
Cancer Institute, 106(5), dju049.
Wang, W. (2003). Estimating the association parameter for copula models under dependent censoring,
Journal of the Royal Statistical Society, Series B, 65(1), 257-273.

コピュラを用いた生存時間解析

173

Wang, Y. C., Emura, T., Fan, T. H., Lo, S. and Wilke, R. (2020). Likelihood-based inference for afrailty-copula model based on competing risks failure time data, Quality and Reliability Engineering International, DOI: 10.1002/qre.2650.
Wu, B. H., Emura, T. and Michimae, H. (2020). Meta-analysis of individual patient data with semicompeting risks under the Weibull joint frailty-copula model, Computational Statistics, DOI:
10.1007/s00180-020-00977-1.
Zhou, Y., Lu, Z., Shi, Y. and Cheng, K. (2019). The copula-based method for statistical analysis ofstep-stress accelerated life test with dependent competing failure modes, Proceedings of theInstitution of Mechanical Engineers, Part O: Journal of Risk and Reliability, 233(3), 401-418.

174

Proceedings of the Institute of Statistical Mathematics Vol. 68, No. 1, 147–174 (2020)

Survival Analysis Using Copulas—Meta-analysis with Correlated Endpoints—
Takeshi Emura1 and Hirofumi Michimae2
1 Department of Information Management, Chang Gung University
2 School of Pharmacy, Kitasato University

With the rapid development this decade of open databases for biomedical researchers,
we have gained access to complex and yet accurate patient-level information. For cancer
patients in particular, these databases record information on individual patients, including overall survival time, time-to-tumor progression, tumor size, and gene expressions.
In addition, there has been signiﬁcant development of meta-analytical methodologies foranalyzing data from diﬀerent sources. To fully utilize such complex survival data, it is insuﬃcient to apply the classical tools used in survival data analysis, such as Cox regression.
In this paper, we review copula-based statistical methods for analyzing two survival time
variables, namely overall survival time and time-to-tumor progression. We also review thejoint frailty-copula model for individual patient data (IPD) meta-analysis methods, which
account for the heterogeneity of patients from diﬀerent sources. Regarding the construction of an appropriate likelihood function for a given dataset, we explain the importanceof acknowledging the semi-competing risks relationship between overall survival time andtime-to-tumor progression. Finally, we introduce a dynamic prediction method for overall
survival time according to gene expressions and tumor progression, which may contributeto the development of personalized medicine.

Key words: Clinical trial, competing risk, Cox proportional hazards model, dynamic prediction, gene expression, personalized medicine.