統計数理（2009）
第 57 巻 第 1 号 119–138
c 2009 統計数理研究所


特集「確率過程の統計解析」［総合報告］

生存時間解析におけるセミパラメトリック推測と
その周辺
服部 聡

†

（受付 2008 年 6 月 27 日；改訂 2008 年 8 月 21 日）

要旨

加法ハザードモデル，比例オッズモデル，加速モデルなど，比例ハザード性を満たさないセ
ミパラメトリックモデルに対する様々な推測法が，近年提案されてきた．Cox 比例ハザードモ
デルで用いられる部分尤度最大化は，これらのモデルには有効に働かないが，多くのモデルに対してマルチンゲール推定方程式の立場から，推測法が統一的に構成可能である．また，マルチンゲール残差に基づいて統一的に回帰診断を行うことができる．本論文では，これらの方法を中心として，生存時間解析におけるセミパラメトリック推測および関連する話題についての最近の進展を概観する．
キーワード： Cox 比例ハザードモデル，加速モデル，加法ハザードモデル，経験過程，
線形変換モデル，マルチンゲール．
1.

はじめに

生存時間解析とは，何らかの関心のある事象（イベント）が生起するまでの時間に関する統計的推測を行う分野の総称であり，医学・経済学・品質管理などの分野で広く応用されている．医学分野については，例えば臨床腫瘍学の専門誌である Journal of Clinical Oncology 誌などを見ると，多くの論文において，生存時間解析分野の基本的な道具である Kaplan-Meier 法，logrank
検定，Cox 比例ハザードモデルが頻繁に用いられている．腫瘍領域の臨床研究では，死亡までの期間，病勢進行までの期間など，状態の悪化までの期間を問題にすることが多く，これらの方法は不可欠な道具となっている．腫瘍領域以外でも，治癒するまでの期間など，時間を問題にすることは多く，生存時間解析が広く応用されている．
このような生存時間（死亡をイベントとしない場合も，このように呼んでおく）
を評価する臨床研究においては，有害事象の発現などにより，イベントが観察される以前に研究から脱落す
ることが不可避である．また，全ての症例のイベントが観察されるまで研究を継続すると，研究結果を得るのに多大な時間を要することから，全ての症例のイベントの観察を待たずに統計解析を行うことが通常である．その場合，一部の症例に対してはイベントが観察されていないことになる．このような症例は（右側）打ち切りと呼ばれ，Kaplan-Meier 法，logrank 検定，Cox
比例ハザードモデルを始めとする生存時間解析分野の諸手法は，何らかの方法で打ち切り症例を考慮して構成される．打ち切りを考慮して構成した諸手法は一般に複雑になり，その理論的考察も標準的なものではなくなる．Kaplan-Meier 法，logrank 検定，Cox 比例ハザードモデル†

久留米大学 バイオ統計センター：〒830–0011 福岡県久留米市旭町 67

統計数理

120

第 57 巻

第 1 号 2009の正当化には，マルチンゲールなどの確率過程論の技術が本質的な役割を果たしており，実用上極めて重要な諸手法に確率過程論がエレガントに応用されている点が，この分野を興味深いものとしている（Fleming and Harrington, 1991; Andersen et al., 1993; Kalbﬂeisch and Prentice,．
2002）
Cox（1972）
により提案された Cox 比例ハザードモデルは，基準ハザード関数をパラメトライ
ズしないセミパラメトリックモデルであり，一般にパラメトリックモデルよりも柔軟であると考えられる．そのため，SAS などの統計ソフトウェアーで手軽に実施できることも相まって，極めて広範に用いられている．一方で，Cox 比例ハザードモデルが要求する比例ハザード性の仮定は，実際のデータ解析の場面では，必ずしも満たされない．その場合，他のモデルを検討することになるが，統計ソフトウェアーで手軽に実行できることから，パラメトリック加速モ
デルが用いられることも多い．しかしながら，その際には基準ハザードをどのようにパラメトライズするかが問題となる（Kalbﬂeisch and Prentice, 2002）．1980 年代後半から，基準ハザード
関数を特定することなく，かつ比例ハザード性を満たさないセミパラメトリックモデルの推測理論が著しく発達してきた．本論文では，Cox 比例ハザードモデルおよび比例ハザード性を要求しないセミパラメトリックモデルの推測理論を中心に，最近の生存時間解析の発展の概観を行う．生存時間解析の推測法の発展には，確率過程論的な考察が不可欠であるが，背景にある
確率過程論的な技術的側面にも簡単に触れることとする．本論文の構成は以下の通りである．
2 節で，本論文が主に取り上げる一変量生存時間解析における基本的な記法を準備する． 3 節
では Cox 比例ハザードモデルの推測について簡単にまとめる．Cox 比例ハザードモデルのパラメータ推定は部分尤度の最大化に基づくが，マルチンゲール推定方程式としての導出も可能であり，基本的な考え方をまとめる．4 節では，比例ハザード性を要求しないいくつかのセミパ
ラメトリックモデルの推測について紹介する． Cox 比例ハザードモデルで用いられた部分尤度の方法は，これらのモデルには有効には働かず，モデルの特質を利用した様々な推測法が提案されている．一方で，多くの場合に，マルチンゲール推定方程式の観点から統一的に推測法を導出でき，その点を中心に解説する．5 節ではこれらのセミパラメトリックモデルに対する回
帰診断技術についてまとめる．特に，多くのセミパラメトリックモデルに対して統一的に回帰診断を行うことができる累積マルチンゲール残差による方法を中心に解説する．6 節ではいく
つかの関連する話題として，縮小回帰モデル，多状態モデル，経時測定データに対する生存時
間解析の応用についての最近の研究について紹介する．まとめと，言及できなかったいくつかの話題について，7 節で言及する．
2.

生存時間データ

ある時点を原点とし，関心のあるイベントが生起するまでの時間を表す確率変数を T とし，打ち切りを示す確率変数を C とする．X = min(T, C) とし，イベントが観察されたか否かを示
す二値確率変数を ∆ = I(T ≤ C) とする．Z を p 次元共変量とし，有界であると仮定する．実際には，多くの場合に時間依存共変量の場合も取り扱うことができるが，統一的に議論を進めるた
めに，時間 t に依存しないもののみを考える． T あるいは C のいずれか一方しか観測できない
右側打ち切りの場合を考え，(X, ∆, Z) が観測されるものとする． Z を与えたもとで，打ち切り
C はイベント T と独立であるとする． T は連続型の非負の確率変数とし，Z を与えたときの Tの条件付き生存時間関数を S(t|Z) = P (T > t|Z)，条件付きハザード関数を λ(t|Z) = −dS(t|Z)/dt，
S(t|Z)
t
累積条件付きハザード関数を Λ(t|Z) = 0 λ(u|Z)du により定義する．T, C などの確率変数に対し，i 番目の症例に対する観測値を Ti , Ci などと下付き数字で表すこととする． i 番目の症例に対し，観測値 (Xi , ∆i , Zi ) が得られているものとし，{(Xi , ∆i , Zi )}ni=1 は独立に (X, ∆, Z) と同

生存時間解析におけるセミパラメトリック推測とその周辺

121

一の分布に従うものとする．本稿ではこの右側打ち切りを伴う一変量生存時間解析の枠組みを基本として議論を進めることとする．その自然な拡張と考えられる再発事象データについても扱うが，必要な記法は後に準備する．
3.

Cox 比例ハザードモデルと部分尤度による推測

3.1 Cox 比例ハザードモデル
Cox（1972）
により提案された Cox 比例ハザードモデルはλ(t|Zi ) = λ0 (t) exp(β0T Zi )により定義される．ここで，λ0 (t) は基準ハザード関数と呼ばれ，β0 は p 次元の回帰係数である．イベントが観察された時間を t(1) < t(2) < · · · < t(L) とする．同一時刻に複数の症例にイベントが生じていないとし，t(l) にてイベントが観察された症例の共変量を Z(l) と書く．R(t) を
時間 t におけるリスク集合とする．つまり，時間 t の直前までにイベントが観察されず，かつ打ち切りも受けていない症例全体とする．β0 の推定は，部分尤度関数(3.1)l(β) =

L
i=1

exp (β T Z(i) )Tl∈R(t(i) ) exp (β Zl )の最大化により行われることが多い（Cox, 1972）．対応するスコア関数は



L
n



T
T
(3.2)
UP H (β) =log
exp (β Zl )β Z(i) −i=1

i=1l∈R(t(i) )により与えられる．UP H (β) = 0 の解を β̂P H とする．部分尤度に基づくスコア関数は基準ハザード関数 λ0 (t) を含んでおらず，基準ハザード関数 λ0 (t) を特定することなく， β̂P H を標準的な Newton-Raphson 法により得ることができる．一方で，スコア関数は独立な確率変数の和の形になっておらず，扱いが面倒なものとなっている．Tsiatis（1981）は独立な確率変数の和に近似することにより，一致性，漸近正規性などの β̂P H の理論的正当性を与えた．その後
Andersen and Gill（1982）が，より一般的な設定のもとで，見通しのよい理論的正当化を与えた．τ を P (Ci > τ ) > 0 を満たす定数とし， Ni (t) = I(Xi ≤ t, ∆i = 1)，Yi (t) = I(Xi ≥ t) とする．
Ni (t) および Yi (t) はそれぞれ計数過程，リスク集合過程と呼ばれるもので，データ (Xi , ∆i ) と{(Ni (t), Yi (t))} は一対一に対応する．したがって，(Ni (t), Yi (t), Zi ) をデータとみなすことがで
きる．Ft = σ{Ni (u), Yi (u), Zi : u ≤ t} とする．ただし，σ(·) は · で生成される最小の σ 代数とし，その族 F = {Ft : t ∈ [0, τ ]} は ﬁltration と呼ばれる．Ni (t) は F に関して劣マルチンゲールになる．このことから，適当な正則条件のもと，F− マルチンゲール Mi (t) が存在して，(3.3)

Ni (t) −t
0

Yi (u)dΛ(u|Zi ) = Mi (t)と一意的に分解できることが知られている（Fleming and Harrington, 1991 など）．ただし，
t
Λ(t|Zi ) = 0 λ(u|Zi )du とする．この分解は Doob-Meyer 分解と呼ばれ，左辺第二項は，可予測過程と呼ばれるものであり，補正項（compensator）と呼ばれる．可予測過程の正確な定義は省略するが，実用上頻繁に現れる可予測過程は，左連続な sample path を持つ確率過程である．
時点 t において，直前 t− までの事象を条件付けると，左連続性により左辺第二項は定数となり，Mi (t) を誤差項とみることで，Doob-Meyer 分解は「データ − 回帰部分=残差」の構造と解釈され，ハザードに対してモデリングすることが，生存時間データに対する自然な回帰モデルの枠組みであることが理解される．

122

統計数理

第 57 巻

第 1 号 2009

部分尤度に基づくスコア関数は
UP H (β) =

n

i=1

τ
0

{Zi − Z̄(u; β)}dNi (u)と書くことができる．ただし，確率変数 Ni (t) に関する積分は，sample path ごとに LebesgueStiltjes 積分を用いて定義するものとし，S (1) (u; β)S (0) (u; η)nYj (u) exp(β T Zj )S (0) (u; β) = n−1Z̄(u; β) =j=1nZj Yj (u) exp(β T Zj )S (1) (u; β) = n−1
j=1とする．この表現と Doob-Meyer 分解（3.3）
を介して，マルチンゲールに関する収束定理である
Rebolledo のマルチンゲール中心極限定理や Lenglart の不等式
（Andersen et al., 1993 など）を応用することで，Andersen and Gill（1982）
は Cox 比例ハザードモデルにおける β̂P H の一致性，
t
漸近正規性を示した．さらに基準累積ハザード関数 Λ0 (t) = 0 λ0 (u)du の Nelson-Aalen 型の推定量である Breslow 推定量（Breslow, 1972）Λ̂(t) =n
dNi (t)(0) (t; β)nSi=1の漸近正規性を示した．
ところで，Andersen and Gill（1982）での議論では（3.3）の分解が本質的であって， Ni (t) がNi (t) = I(Xi ≥ t, ∆i = 1) により定義されていることは本質的ではない．一人の症例から複数のイベントが繰り返し観察されうる場合を想定し，Ni (t) を時間 t までに症例 i に生じた総イベン
ト数，つまりイベントが生じるごとに 1 だけ跳躍する確率過程と考え，Yi (t) を症例 i が時間 t
においてリスクにあれば 1，無ければ 0 をとる確率過程とし，ハザード関数を強度関数
λ(t|Zi ) = lim h−1 P (Ni (t + h) − Ni (t) = 1|Ft )h↓0に取り替えることで，（3.3）を通じて再発事象に関する Cox 型の回帰モデルを考えることができる．Andersen and Gill（1982）はこの枠組みで一気に Cox 回帰モデルに対する部分尤度による推測法の正当化を与えており，再発事象に対する彼らのこの枠組みは Andersen-Gill モデルと呼ばれることがある．彼らのマルチンゲールによる理論的正当化は，ハザード関数あるいは強度関数が Cox 回帰型であることに強くは依存しておらず，次節で議論するセミパラメトリッ
クモデルなど，様々な状況に適用可能である．もともと，マルチンゲールによる方法は，Aalen（1978）により一標本問題における累積ハザード関数に対する Nelson-Aalen 型推定量の性質を調べるために用いられ，Gill（1980）
は，この方法により，二標本問題における logrank 検定，一般化 Wilcoxon 検定を含む検定統計量の性質を調べた．このように，マルチンゲールの方法は生存時間解析の極めて広範な問題に対して統一的に適用可能であり，生存時間解析研究の日常の道具として広く用いられている．
3.2 部分尤度による推測のマルチンゲール推定関数としての解釈
部分尤度（3.1）は，全尤度を適当な条件付き確率密度の積に分解し， β0 に関する情報をほとんど含まない部分を無視することにより得られ（Cox, 1975; Fleming and Harrington, 1991），そ生存時間解析におけるセミパラメトリック推測とその周辺

123の結果 Cox 比例ハザードモデルの場合には，基準ハザード関数に依存しないスコア関数（3.2）
を得ることができる．このように，Cox 比例ハザードモデルの推測は尤度によるものと考えられるが，異なる立場からスコア関数（3.2）を導入することができる．
Mi (t; β, Λ) = Ni (t) −t0

Yi (u) exp (β T Zi )dΛ(u)とし，次のようなマルチンゲールに基づく推定方程式のペアを考える．
n

Mi (t; β, Λ) = 0

i=1
n

i=1

τ
0

Zi dMi (u; β, Λ) = 0.

E[Mi (t; β0 , Λ0 )|Zi ] = 0 により，この推定方程式を考えることは自然である．第一式が全ての tに成り立つことから，n
i=1

dNi (u) =nYi (u) exp(β T Zi )dΛ(u)

i=1

n

i=1 dNi (u)TY
i=1 i (u) exp(β Zi )⇐⇒ dΛ(u) = nとなる．これを第二式に代入して，dΛ(u) を消去するとスコア関数（3.2）
が得られる． β を β̂P Hとしたとき，Λ(u) は Breslow 推定量に一致する．したがって，この導出は，初めに部分尤度により β̂P H を得て，その後 Breslow 推定量を得る通常の流れの逆を行っていることになる．本
論文では，この方法をマルチンゲール推定方程式による方法と呼ぶこととする．次節以降で見るように，この方法は広範なセミパラメトリック推測に適用可能な方法となっている．
3.3 平均値関数に対する Cox 型のモデルへの拡張
Andersen-Gill モデルは Doob-Meyer 分解の成立が前提となっており，それは {Ni (t)} がポア
ソン過程であることを要請することに相当する． Lin et al.（2000）は，この仮定を緩め， Ni (t)の平均値関数に Cox 型のモデルを仮定したモデル
E[Ni (t)|Zi ] = Λ0 (t) exp(β0T Zi )に対して，通常の部分尤度に基づくスコア関数による推定量が，一致性および漸近正規性を有することを示した．この枠組みでは，Doob-Meyer 分解が成立せず，従ってマルチンゲールによる理論的な議論は利用することができないことになる．Lin et al.（2000）は経験過程の理論（Pollard, 1990; van der Vaart and Wellner, 1996; Bilias et al., 1997）の応用により，理論的正当化を行った．
4.

比例ハザード性を要求しないセミパラメトリックモデルの推測

生存時間解析に Cox 比例ハザードモデルが与えた影響は大きく，基準ハザード関数を特定しないセミパラメトリック推測がデータ解析法の中心的な役割を果たしている．一方で，比例ハザード性を始めとするいくつかの強い仮定を Cox 比例ハザードモデルは要求しており，実際のデータ解析において，必ずしも適切ではない．そのため，基準ハザードを特定せず，比例ハ
ザード性を要求しないセミパラメトリックモデルの推測が発展してきた．本節では，そのようなセミパラメトリックモデルを概観する．モデルの定式化を巧妙に利用した方法が多く提案されている一方で，統一的な視点で推測法を理解することも可能となっている．

統計数理

124

第 57 巻

第 1 号 2009

4.1 加法ハザードモデル
加法ハザードモデルはλ(t|Zi ) = λ0 (t) + β0T Zi(4.1)

により定義される（Breslow and Day, 1980, p. 183）
．このモデルに対して，部分尤度による推測を考えても，基準ハザード関数に依存しない β0 の推定方程式を導入することは困難である．Linは以下のような基準ハザード関数に依存しない β0 の推定方程式を提案した．
and Ying（1994）(4.2)

UAH (β) =ただし，Z̄(t) =n
j=1 Yj (t)Zjn
j=1 Yj (t)の解 β̂AH は(4.3)

β̂AH =

 n

i=1

τ
0

n

i=1

τ
0

Zi − Z̄(t)dNi (t) − Yi (t)β T Zi dtとする．この推定方程式は明示的に解くことができ， UAH (β) = 0

Yi (t) Zi − Z̄(t)⊗2dt−1  n

i=1

τ
0


Zi − Z̄(t) dNi (t)により与えられる．
は部分尤度に基づくスコア関数としての解釈はされないが，マル
ところで，推定方程式（4.2）
チンゲール推定方程式として導くことができる．
t

Mi (t; β, Λ) = Ni (t) −0

Yi (u){dΛ0 (u) + β T Zi du}として推定方程式のペアn

Mi (t; β, Λ) = 0

i=1
n

i=1

τ
0

Zi dMi (u; β, Λ) = 0

を考える．第一の式より
n

dNi (u) =

i=1

n


Yi (u)dΛ(u) + β T

i=1⇐⇒ dΛ(u) =

n

n


Yi (u)Zi du

i=1

i=1 dNi (u) − Yi (u)βni=1 Yi (u)

T

Zi du

とし，これから第二式の dΛ(u) を消去することで，（4.2）を得ることができる．
打ち切りがない場合には，線形回帰モデルは明示的な解を有しているが，後述するように，
線形回帰モデルに対応する加速モデルは，打ち切りが存在する場合には，推定量は明示的な解を持たない．一方で，加法ハザードモデルの定義（4.1）の右辺は負の値をとり得るという欠点を有しているが，（4.3）のように推定量が明示的に与えられる利点を有している．本論文では詳し
くは触れないが，区間打ち切りと呼ばれる打ち切りを受ける生存時間データなど複雑なデータに対して，簡明なセミパラメトリック推測法を構成する方法として，加法ハザードモデルが有効な方法となっている（Lin and Ying, 1997; Martinussen and Scheike, 2002a）．
4.2 比例オッズモデル・線形変換モデル
比例オッズモデルはlog

1 − S(t|Zi )1 − S0 (t)= log+ β0T Zi
S(t|Zi )S0 (t)

生存時間解析におけるセミパラメトリック推測とその周辺

125にて定義される（Bennett, 1982）．ただし S0 (t) は基準生存時間分布とする．回帰係数 β0 の各成分は，説明変数の違いにより， t 年生存率に対するオッズに対して t に依存せずに定数倍される効果を示している．この比例オッズ性は，t → ∞ によりハザード比が 1 に収束することが示され，ハザード比の意味で説明変数の影響が減弱していく状況を反映したモデルということになる．比例オッズモデルに対する推測としては，Yang and Prentice（1999），Murphy et al.
（1997），Scharfstein et al.（1998）などがある．
線形変換モデルはg(S(t|Zi )) = h0 (t) + β0T Zi(4.4)

により定義される．ただし，h0 (·) は未知の単調非減少関数であり，g(·) は既知の単調非増大関数である．g(x) = log 1−xとした場合が，比例オッズモデルであり， g(x) = log (− log x) とした
x
場合が，Cox 比例ハザードモデルに相当する．
（4.4）において t = T とすることで，線形変換モ
デルはh0 (Ti ) = −β0T Zi + i

(4.5)

と表現することができる．ただし，i は既知の分布に従う確率変数で，i が極値分布に従う場
合が，Cox 比例ハザードモデルに相当し，標準 logistic 分布に従う場合が，比例オッズモデルに相当する．この表現を巧妙に利用し，Cheng et al.（1995），Fine et al.（1998）
はノンパラメト
リック部分 h0 (t) に依存しない β0 に対する U 統計量型の推定方程式を提案し，一致性・漸近正規性を示した．
Chen et al.（2002）
は，マルチンゲール推定方程式の考え方による推測法を提案した．推定方
程式の組
n


[dNi (t) − Yi (t)dΛ{h(t) + β T Zi }] = 0

i=1
n

i=1

τ
0

Zi [dNi (t) − Yi (t)dΛ{h(t) + β T Zi }] = 0を考える．Cox 比例ハザードモデル，加法ハザードモデルにおいては，固定した β に対して第一の方程式を基準ハザードに関して明示的な式により解いていたが，線形変換モデルでは明示的に解くことはできないことになる． Chen et al.（2002）
は，この二つの推定方程式を，基準ハザードおよび β について交互に解くことを提案し，その結果得られる推定量が，一致性および漸近正規性を持つことを示した．
4.3 加速モデル
加速モデルは(4.6)

log Ti = β0T Zi + i

(4.7)⇐⇒ Ti = exp (i ) exp (β0T Zi )により定義されるモデルである．ここで，i は Zi と独立な確率変数で，その確率分布関数 F0 (·)
は未知であるものとする．このモデルは（4.7）より，基準ハザードに対応する確率分布 F0 (·) にT従う確率変数 exp (i ) を exp (β0 Zi ) exp (i ) に加速させていると解釈される．加速モデルの条件
付ハザード関数は
Λ(t|Zi ) = Λ0 (t exp(−β0T Zi ))

126

統計数理

第 57 巻

第 1 号 2009を満たす．ただし Λ0 (t) = − log{(1 − F0 (t))} は，基準ハザード関数とする．すなわち，加速モ
デルは説明変数により時間軸のスケール変換を行っていることに相当する．
定義（4.6）により，加速モデルは誤差項の分布を特定しない線形回帰モデルであるので，打ち切りが生じない場合には，最小二乗法により明示的な推定量が与えられる．しかしながら，打ち切りが存在する場合には，明示的な解を得ることができず状況は一変する．Buckley andは，右側打ち切りが存在する場合の最小二乗推定量である Buckley-James 推定量
James（1979）を提案した．彼らの方法は，打ち切り症例のイベントを，データを与えたもとでの条件付期待値で代用し，最小二乗法を適用する方法であり，再帰的な計算が必要となる．更にその理論的性質を導出することも容易ではなくなる．Lai and Ying（1991）は，chaining などの確率過程論
的な技術を用いて一致性・漸近正規性を示した． Jin et al.（2006）
は次に述べる順位回帰による方法に対して Jin et al.（2003）が提案した線形計画法および標本再抽出法による比較的簡明な方法の提案を行った．
次に，Prentice（1978）に始まる順位回帰に基づく方法について述べよう． i (β) = log Xi − β T Ziと定義する．{i (β), ∆i } に対する logrank 検定統計量1
∆i ψ(β, ei (β))(Zi − Z̄i (β))
n i=1
n

(4.8)

UAF T (β) =を考える．ただし，ψ は重みとし（Jin et al., 2003），
nj=1 Zj I(ej (β) ≥ ei (β))
Z̄i (β) = n
j=1 I(ej (β) ≥ ei (β))
ei (β) = log(Xi ) − β T Ziとする．UAF T (β0 ) の各成分は一標本での（重み付き）logrank 検定統計量であることから，漸近的に平均 0 の正規分布に従う．このことから，逆に logrank 検定統計量が 0 となる β により，
β0 が推定されることが期待される．実際には（4.8）
は階段関数であり，一般にはゼロ点を有さないので， UAF T (β)  を最小にする β̂AF T により推定される．  UAF T (β)  は β に関して滑らかでないことから，その最小化はやはり標準的ではないが，Lin and Geyer（1992）は，simulated
annealing による方法を開発した．また，推定方程式が滑らかでないことから，漸近正規性の証
明は標準的なものではなくなるが，Tsiatis（1990）, Lai and Ying（1992）, Ying（1993）などによりなされた．漸近分散は未知の誤差項 i の確率密度関数 dF0 (t)/dt に依存しており，直接漸近分散を推定するには平滑化を用いる必要がある． Wei et al.（1990）は， minimum dispersion 統
計量により，確率密度関数 dF0 (t)/dt を推定することなしに，回帰係数の検定を行う方法を提案した．回帰係数の次元が小さい場合には，この検定統計量に基づいて信頼区間を構成することが可能であるが，一般には困難となる．Lin et al.（1998）, Jin et al.（2001）は，推定方程式を標本再抽出することで，確率密度関数 dF0 (t)/dt を推定することなしに，回帰係数の信頼区間を構成する方法を開発した．その標本再抽出の方法は， Random Multiplier 法と呼ばれ，後述
する累積残差に基づく回帰診断，生存時間関数の同時信頼区間の構成（Lin et al., 1994; Cheng，分位点回帰（Parzen et al., 1994; Jin et al., 2001）などの，幅広い問題で応用されet al., 1997）ている．一般に β̂AF T は一意ではないが，基礎となる logrank 検定統計量として Gehan スコア（Gehan, 1965）
を用いると，推定方程式は凸関数となり，漸近的にを与えた logrank 検定統計量
一意的な解を有する．Jin et al.（2003）
は，重みに Gehan スコアを与えた logrank 推定方程式による推定値が線形計画問題を解くことで得られることを示し， Random Multiplier 法により標本再抽出した推定方程式を線形計画法により解くことで効率的に信頼区間を構成する方法を提案した．更に，広範なスコアに対して，Gehan スコアによる推定値を初期推定値として再帰的

生存時間解析におけるセミパラメトリック推測とその周辺127に推定する方法も提案した．
このように，基本的とも考えられる線形回帰問題ではあるが，打ち切りがある場合には困難が多く，非常に豊かな統計的方法が発展している領域である．最近，セミパラメトリックモデ
ルの意味での漸近有効な推定法が Zeng and Lin（2007a）
により提案された．logrank 推定方程式
による方法は自然に Andersen-Gill 型の加速モデルに拡張することができ（Lin et al., 1998; Jin，打ち切りがある場合の線形回帰モデルの推測として順位回帰による方法は自然な
et al., 2005）
方法であるといえる．
加速モデルは，定義（4.7）
より，生存時間データを，共変量 Zi に依存しない，従って各症例に共通な分布を有する確率変数に変換するモデルと解釈される．このように，モデルの定義が確率変数に対して直接的に行われているのは，Cox 比例ハザードモデルなどの他のモデルにはない性質であり，この確率変数間の変換を，因果推論における反事実仮想下への変換と見なすことで，無作為化臨床試験における治療不遵守の調整において，加速モデルは重要な役割を果た
している（Robins and Tsiatis, 1991; Robins, 1993）．また，繰り返し測定データの解析において，
いわゆる情報を持つ脱落（Diggle et al., 2002 など）が生ずる場合には，推定にバイアスが生じるが，この状況においてバイアスのない推定を可能にする方法のひとつである artiﬁcial censoring法においても加速モデルは重要な役割を果たしている．この方法は，情報を持つ脱落時刻を加速モデルにより変換し，変換された脱落時刻以降のデータを人工的に脱落させる方法であり，加速モデルの確率変数間の変換という解釈が重要な役割を果たしている（Lin and Ying, 2003）．
5.

セミパラメトリックモデルに対する回帰診断法

5.1 適合度検定
Cox 比例ハザードモデルに対しては様々な観点からの適合度検定の方法が提案されている．
Lin and Wei（1991）
は，Cox 比例ハザードモデルによる推定量の分散が，部分尤度を尤度関数と見なしたときの Fisher 情報量の逆数で与えられることを利用して，二種類の情報行列の差を評価する White（1982）による情報検定に対応する適合度検定を提案した． Jones and Harrington（2001）
はマルチンゲールが平均 0 でかつ独立増分であるという性質に着目した適合度検定を提案した．Lin（1991）は，通常の部分尤度による推定量と，部分尤度に基づくスコア関数に重みを導入した，重み付きスコア関数による推定量が，Cox 比例ハザードモデルを正しく特定した場合には同一の真のパラメータ β0 に収束するものの，誤特定した場合には異なる極限を有することを利用した適合度検定の方法を提案した．この方法は Gill and Schumacher（1987）
により，二標本の場合の比例ハザード性の適合度検定として導入されたが，重みつきの推定方程式を考えることができる状況であれば適用でき，実際加速モデルの適合度検定にも応用された（Wei et．これらの方法は理論的にも妥当な適合度検定であり，モデルの当て
al., 1990; Lin et al., 1998）
はまりを客観的に評価することができるが，一方で，このような包括的な適合度検定に共通の問題として，一旦モデルの適切性が疑われた際に，どのような誤特定が生じているかという，モデルを改善するための情報は与えられないことが挙げられる．
5.2 累積残差による回帰診断法
回帰モデルの診断において，残差を調べることは非常に有用であるが，生存時間解析における回帰モデルに対しては複数の残差の提案がある．定義が自然であり，広く用いられている残差として，マルチンゲール残差がある． Cox 比例ハザードモデルの場合，マルチンゲール残差は(5.1)

M̂i (t) = Ni (t) −t0

Yi (u) exp (β̂ T Zi )dΛ̂(u)なる確率過程として定義される．ただし，通常マルチンゲール残差と言った場合には，この定

128

統計数理

第 57 巻

第 1 号 2009

義で t = ∞ としたM̂i = ∆i − Λ̂(Xi ) exp(β̂ T Zi )のことを指す場合がある（Therneau and Grambsch, 1990）．実際統計解析ソフトウェアー SASの PHREG プロシージャではこれが出力される．本稿では（5.1）をマルチンゲール残差と呼ぶ．
Doob-Meyer 分解（3.3）を考えると，この残差がもっとも自然なものであると考えられる．
回帰診断において，残差をプロットすることは回帰モデルに含めた説明変数の関数形の適切性を確認したり，外れ値が存在するか否かを検討する際に有用となる．実際，正規分布の誤差を仮定した重回帰モデルなどでは，残差を各説明変数方向にプロットし，系統的なパターンが見出されないか，あるいは著しく外れた値がないかなどを検討することがよく行われる．このような方法が有効であるには，残差がどのようにばらつくかに関する知識が必要で，誤差の正規性の仮定が重要となる．しかしながら，マルチンゲール残差は，値域が (−∞, 1] と著しく
0 を中心に非対称となっており，単純にマルチンゲール残差をプロットしても，そこからモデ
ルの適切性を判断することは非常に困難である．モデルが正しく特定されていれば，マルチン
ゲール残差は，Doob-Meyer 分解におけるマルチンゲールを近似するものであることから，平均は 0 であることが期待できる．このことから Therneau and Grambsch（1990）, Grambsch et al.
（1995）
はマルチンゲール残差のプロットを平滑化する方法を提案した．しかしながら，結果が平滑化の方法に依存しているうえ，p 値などによるばらつきを考慮した評価ができないことから，主観的であるという欠点を有する．これらの欠点を克服する方法として，Lin et al.（1993）
は累積残差による回帰診断法を提案した．検定統計量として
(5.2)

1

H(t, z) = n− 2nI(Zi ≤ z)M̂i (t)i=1を考える．ただし，z は p 次元ベクトルであり，I(Zi ≤ z) は Zi の各成分がそれぞれ z の各成分よりも小さいときに 1，いずれかの成分で成り立たないときに 0 をとるものとする．モデル
が正しく特定されているとき，この多次元パラメータで添え字付けられた確率過程が，平均 0の正規確率過程に弱収束することが示される．従って，観測された sample path が平均 0 の正規確率過程の sample path でない場合には，モデルの誤特定が疑われることになる．そのためには，極限の正規確率過程の共分散構造を知る必要があるが，解析的にその構造を求めることは困難である．Lin et al.（1993）は random multiplier 法と呼ばれる標本再抽出法により，乱数を発生させることにより，この正規確率過程の sample path を計算機上で発生させる方法を提案した．これにより，モデルが正しいという仮説のもとでの sample path を任意の個数発生させることができ，Kolmogorov-Smirnov 型の統計量supt∈[0,τ ],z∈[0,1]|H(t, z)|の p 値が計算でき，フォーマルな検定が可能となる．更に，
1

H (k) (z (k) ) = n− 2

n


(k)

I(Zi

≤ z (k) )M̂i (τ )

i=1

や
1

H exp (x) = n− 2nI(β̂ T Zi ≤ x)M̂i (τ )

i=1
(k)

といった統計量を考えることができる．ただし，Zi を Zi の第 k 成分とし，z (k) を z の第 k成分とする．これらも，当てはめたモデルが正しいときに，平均 0 の正規確率過程に弱収束す生存時間解析におけるセミパラメトリック推測とその周辺129ることが示され，やはり random multiplier 法により帰無分布を標本再抽出法により発生させることができる．前者は k 番目の説明変数に沿って残差を集めていることに相当し，したがって，k 番目の説明変数の関数形の誤特定に対する指向性の強い検定となる．また後者は比例ハ
ザードモデルにおいて， β0T Zi を指数関数での変換によりモデル化している部分の仮定の誤特定に対して指向性の強い検定と考えられる．これらの統計量は一次元の確率過程であることか
ら，観測された sample path を random multiplier 法により計算機上で発生させた帰無仮説下での sample path と共にグラフィカルに表示することで，両者の違いを視覚的に検討することができる．これは，標準的な重回帰モデルにおいて，残差を特定の説明変数方向にプロットしていることに相当しており，更にマルチンゲール残差の非対称性による解釈の困難さを，random
multiplier 法による帰無仮説下での sample path とグラフィカルに比べることで回避していることになる．このような様々な指向性を持つ検定統計量を検討することで，モデルの誤特定が示唆される場合に，どのような誤特定が生じているかのヒントを得ることができ，非常に有用である．
Lin et al.（1993）
は一変量の Cox 比例ハザードモデルに対してこの方法を提案したが，平均値関数に対する Cox 型のモデルに対しても適用可能である（Lin et al., 2000）．また，マルチ
ンゲール残差は自然に他のモデルにも定義され，累積残差による回帰診断が適用可能となる．
Yuen and Burke（1997），Kim et al.（1998）は加法ハザードモデルに対して，同様の方法を提案し，Hattori（2008a）は，線形変換モデルに対しての累積残差による方法の正当性を示した．
加速モデルに対しては，Lin and Spiekerman（1996）はパラメトリック加速モデル（Kalbﬂeischの回帰診断法を累積マルチンゲール残差により提案した．Stute et al.
and Prentice, 2002 など）（1998）, Stute（2000）は，セミパラメトリック加速モデルに対して，累積残差による方法を提案した．ただし，残差としてはマルチンゲール残差ではなく，定義（4.6）における i を用いており，Wild bootstrap 法と呼ばれる標本再抽出法を用いている．累積残差による方法は生存時間解析に特有の方法ではなく，一般化線形モデル（Su and Wei, 1991）にも適用でき，繰り返し（Pan and Lin,
測定データに対する一般化推定方程式（Lin et al., 2002），一般化線形混合モデル，セミパラメトリックモデル（Hattori, 2008b）など，広範な問題に対して適用されている．
2005）
累積残差による方法により， Kolmogorov-Smirnov 型の検定統計量に基づいて，フォーマルな検定が可能となり，更に，指向性のある検定統計量と帰無分布をグラフィカルに表示することで，検定統計量の sample path が帰無分布のそれと比べて奇妙か否かの視覚的な判断ができることから，非常に有意義ではある．一方で，そのグラフ表示の結果の解釈は，残差を累積していることから直接的ではなく，必ずしも容易ではない．Lin et al.（2002）
は一般化推定方程式の場合ではあるが，様々な誤特定に対して，累積残差がどのような挙動をするかのプロトタイ
プをシミュレーションにより求め，そのプロトタイプを参考にどのような誤特定が生じているかのヒントを得る方法を検討しているが，十分に実用的なものにはなっていないように思われる．更には，残差のプロットによる外れ値の確認も，累積残差による方法では困難である．
6.

いくつかの関連する話題

4 節でまとめたように，セミパラメトリック回帰モデルの推測理論はすでに多くの提案がなされている．しかしながら，通常の回帰分析の単純な適用では意味のある結論が引き出せない状況が存在する．本節ではそのような例として，縮小ランク回帰および多状態モデルについて，
簡単に概観する．生存時間解析の興味深い応用として，繰り返し測定データに対する計数過程
アプローチについても，簡単に紹介する．これらの問題に対しても，これまでに紹介したセミ
パラメトリック回帰モデルならびにマルチンゲール推定方程式の考え方が，重要な役割を果た

130

統計数理

第 57 巻

第 1 号 2009

している．
6.1 縮小ランク回帰モデル
マイクロアレイ技術により得られる遺伝子発現量データ，あるいは EGFR，Her2，YB-1 な
どの分子標的マーカーなどのように，その間に pathway が介在し，関連する変数と生存期間の関連を明らかにしたい状況にしばしば遭遇する（Fujii et al., 2008）．このような場合，関連する変数を同時に含めた回帰モデルによる解析は，他の変数を条件付けた上での解析となり，解釈が容易でない．寧ろ，変数間の pathway を同定し，pathway の生存期間への影響を評価する方が自然と考えられる．主成分回帰あるいは部分最小二乗法による方法は，このような pathway
が介在する場合の解析法として重要と考えられる． Cox 回帰モデルによる方法が Nguyen andRocke（2002），Li and Gui（2004），Li and Li（2004）により，加速モデルによる方法が Huang
and Harrington（2005）
により提案された．Ma et al.（2006）
は，Lin and Ying（1994）
による加法
ハザードモデルに対する推定量が，通常の線形回帰のように明示的な解を持つことを利用し，線形回帰の場合に提案されていた主成分回帰の理論を，生存時間解析に拡張した．Wei and Li（2007）は，pathway-based 回帰の名の下に，遺伝子などの間に存在する pathway のハザードへの影響を評価するための Cox 回帰モデルの推測を boosting により行っている．
6.2 多状態モデル
ここまでで議論した生存時間データはイベントの種類が死亡など一つの場合であった．しかしながら，イベントが生じるまでに他の医学的に関心のある事象を生じた後イベントを生ずる症例が存在するなど，各症例の経過は一般に様々である．Andersen et al.（2000）で扱われた，肝硬変症例を対象とした PROVA 試験では，出血を起こし，その後に死亡する症例と，出血を起こさず死亡する症例が見られた．試験開始時の状態を 0，出血を起こした状態を 1，死亡した状態を 2 とし，薬剤ならびに背景要因が，出血を起こすまでの期間に影響するか，出血を伴わず死亡するまでの期間に影響するか，出血した後に死亡するまでの期間に影響するか，に関心があった．このように複数の状態をイベントと見なしたデータに対する関心が最近高まっており，そのためのモデルは多状態（multi state）モデルと呼ばれる．中間に発生するイベントを考慮することで，疾患の進行を考慮した考察が可能となる．肝硬変の例は，3 つの状態で，0 → 1，
1 → 2，0 → 2 の 3 つの推移を許すモデルであり， illness-death モデルと呼ばれる．様々な例がHougaard（1999）にある．
一般に，{0, 1, . . . , k} の k 種類の状態をとり得るとし，i 番目の症例の時刻 t における状態を
Xi (t) で表すこととする．Xi (t) は非斉時的マルコフ過程であるとする．共変量 Zi を与えたと
きの，時刻 s において状態 h である症例が，時刻 t において状態 j に推移する推移確率を
Phj (s, t|Zi ) = P (Xi (t) = j|Xi (s) = h, Zi )により定義する．状態 h から状態 j に推移する条件付き推移強度関数を
αhj (t|Zi ) = lim

1

∆→0 ∆

P (Xi (t + ∆) = j|Xi (t) = h, Zi )により定義し，各条件付き推移強度関数に Cox 型のモデル
T
αhj (t|Zi ) = αhj0 (t) exp(βhjZi )を仮定する．ただし，αhj0 (t) は推移 h → j に対する基準推移強度関数とし，βhj は回帰係数とする．βhj の推定は部分尤度関数の最大化により行うことができる． Nhji (t) を時刻 t までに生じた h → j の推移の回数により定義される計数過程とし，Yhi (t) = I(Xi (t−) = h) とする．部分
尤度を最大にする βhj は，推移 h → j のみをイベントと見なし，Yhi (t) をリスク集合過程とした

生存時間解析におけるセミパラメトリック推測とその周辺

131

場合の Andersen-Gill モデルと見なして推定することに他ならない（Andersen et al., 1993）．したがって，SAS や S-plus などの標準的なソフトウェアーで比較的容易に推定量の計算が可能である（Andersen and Klein, 2002）．Andersen-Gill 型のモデルであれば任意のモデルを用いることができ，Shu and Klein（2005）は加法ハザードモデルによる多状態モデルの推測を議論した．
一方で，いずれのモデルを用いても，各推移が競合的に生じることから，各推移に対する回帰パラメータの解釈は，一般には容易ではない．解釈が容易な量として，共変量 Z0 を有する症例
が時刻 t に状態 j にいる確率，状態占有確率（stage occupation probability）Pr(Xi (t) = j|Z0 ) を考えることができる．各推移に対するモデルから，積積分（product-integral）
により推定されるが（Andersen and Klein, 2002），これは Z0 に関して複雑な非線形関数となり，共変量の効果としては，解釈が困難となる．Scheike and Zhang（2007）, Andersen and Klein（2007）は，推移強度に対してではなく，状態占有確率に対して，一般化線形モデル型のモデルを考えた．これらの方法では，状態占有確率に対する共変量の影響の評価は直接的であり，更にマルコフ性の仮定を必要としない．一方で，多状態モデルのもともとの動機となる病態間の推移は推定されな
いことになる．このようにいずれの立場に立っても，必ずしも解釈が容易でない部分があり，より解釈の容易な方法の開発が必要と考えられる．
6.3 繰り返し測定データの解析
多くの医学研究においては，関心のある変数は経時的に繰り返し評価され，繰り返し測定デー
タと呼ばれる．繰り返し測定データに対しては，様々な解析手法が提案されてきている（Diggleet al., 2002）．最近，Lin and Ying（2001）, Martinussen and Scheike（1999, 2001）, Cheng andなどにより，計数過程に基づく興味深い方法が提案された．ここでは Cheng and Wei
Wei（2000）（2000）で扱われたモデルを考える．Xi (t) を時刻 t における，関心のある変数の測定値とし，モデル
E{Xi (t)|Zi } = µ0 (t) exp(β0T Zi ),i = 1, 2, . . . , nを考える．ただし µ0 (·) は未知関数とし，β0 は p 次元回帰係数とする．二値データに対しては，このモデルは相対リスクに対する回帰モデルに相当する．経時測定データを評価する際，一般には予定された測定時点は遵守できず，測定時期は不規則となる．症例 i の時刻 t までの測定
回数を Ni (t) とする．すなわち，Ni (t) は測定をイベントと見なした際の計数過程とし，Cox 型のモデル E[Ni (t)|Wi ] = Λ0 (t) exp(η0T Wi ) を仮定する．ただし，Wi は Zi のいくつかの成分からなるベクトルとし，η0 の推定は通常の部分尤度に基づく推定量 η̂ により行うことができる．時t
刻 t までの累積測定値 0 Xi (u)dNi (u) を考え，残差MiX (t; β, η) =t
0{Xi (u)dNi (u) − exp (β T Zi ) exp (η T Wi )Yi (u)dΓ̂(u; β, η)}を考える．ただし，
Γ̂(t; β, η) =

n

i=1

t
0

Yi (u)
X (u)dNi (u)T Z } exp {η T W } i
Y
(u)
exp{β
j
j
j=1 j

n

とする．生存時間解析におけるマルチンゲール推定関数の場合と同様にして，β に関する µ0 (t)に依存しない推定方程式を得ることができる（Cheng and Wei, 2000; Hattori, 2008b）. Lin and, Martinussen and Scheike（1999, 2001）, Sun and Wu（2005）
は加法的モデルに対しYing（2001）て，同様な方法により，簡明な推定方程式を導いている．このように，計数過程の方法は，実際の経時測定データに現れるような，評価時期が共変量に依存する状況へのモデリングとして有用である．また，マルチンゲール推定方程式による方法が，経時測定データの解析に対する

統計数理

132

第 57 巻

第 1 号 2009

セミパラメトリック推測に対して，簡明な方法を提供している．理論的正当化には，経験過程の理論やマーク付き点過程の理論が有効に用いられている．
7.

おわりに

本論文では，右側打ち切りを受ける一変量生存時間データならびに再発事象データに対するセミパラメトリック推測法の最近の進展についてまとめた．Cox 比例ハザードモデルに加えて，
加法ハザードモデル，比例オッズモデル，線形変換モデル，加速モデルを取り上げたが，更にそれらを組み合わせたモデル（Lin and Ying, 1995; Chen and Jewell, 2001; Martinussen and Scheike,
や，一部の回帰係数を時間依存にしたモデル（Martinussen et al., 2002; Martinussen and2002a）の推測法も開発されている．これらの多くでマScheike, 2002b; Scheike and Martinussen, 2004）
ルチンゲール推定方程式の考え方が有効に用いられている．クラスター生存時間データあるいは多変量生存時間データに対する周辺モデルに対しても，本論文で述べた方法は自然に拡張されている（Spiekerman and Lin, 1998; Yin and Cai, 2004; Lu, 2005; Jin et al., 2005 など）．また，
本稿では簡単のため共変量 Zi は時間に依存しないものとして話を進めたが，4 節で扱った推定法の多くで，時間依存型共変量を扱うことができる．
マルチンゲールは生存時間解析の手法の構成ならびに理論的発展に重要な役割を果たしてきた．しかしながら，漸近理論の道具として見た場合，マルチンゲール中心極限定理が有効に働くのは，一変量生存時間の場合あるいはポアソン性を仮定した Andersen-Gill モデルの場合に限られ，3 節に示したように，ポアソン性を要請しない場合や，クラスター生存時間データあるいは多変量生存時間データに対してはうまく機能しない．また，5 節で議論した累積マルチ
ンゲール残差による適合度統計量（5.2）は多次元パラメータで添え字付けられており，マルチ
ンゲールによる理論的正当化はできない．これらの問題に対しては経験過程の理論により正当化が行われている．本論文ではマルチンゲール推定関数に焦点を絞り議論したが，その起源は
Cox 比例ハザードモデルに対する部分尤度法に見出すことができる．一方で，Cox（1972）に対する討論において Breslow（1972）はノンパラメトリック最尤法による部分尤度法の解釈を与えた．その考えに則り，最近，ノンパラメトリック最尤法により，混合モデルを含む広範なセミパラメトリックモデルの推測法が提案され，その理論的正当化ならびに推定量を計算するための EM アルゴリズムが実装された（Zeng and Lin, 2007b）．理論的正当化には，やはり経験過程の理論が応用されており，経験過程の理論は生存時間解析研究に不可欠な道具となっている．
6 節で概観したように，分子標的マーカー，遺伝子などの間に存在する pathway に基づく生存時間解析あるいは，複雑な病態推移を考慮した生存時間解析など，より生物学的，医学的な理解を促進するようなモデルの開発および推測法の開発の重要性が近年増していると考えられる（Aalen and Frigessi, 2007）．経験過程などの理論的基盤により，応用上重要な問題に対する推測法の開発の更なる発展が望まれる．
謝辞

査読者ならびに編集委員の先生方には，誤り，不備をご指摘いただき，有益なご意見をいた
だきました．また，久留米大学大学院医学研究科 室谷健太氏は，初稿を通読し，多くの有益な意見をくださいました．ここに謝意を表します．

生存時間解析におけるセミパラメトリック推測とその周辺

133

参 考 文 献Aalen, O.（1978）. Nonparametric inference for a family of counting processes, The Annals of Statistics, 6, 701–726.
Aalen, O. and Frigessi, A.（2007）. What can statistics contribute to a causal understanding?, Scandinavian Journal of Statistics, 34, 155–168.
. Cox’s regression model for counting process: A large sample
Andersen, P. K. and Gill, R. D.（1982）
study, The Annals of Statistics, 10, 1100–1120.
. Multi-state models for event history analysis, Statistical MethAndersen, P. K. and Klein, J.（2002）
ods in Medical Research, 11, 91–115.
. Regression analysis for multistate models based on a pseudoAndersen, P. K. and Klein, J.（2007）
value approach, with applications to bone marrow transplantation studies, Scandinavian Journal of Statistics, 34, 3–16.
Andersen, P. K., Borgan, O., Gill, R. D. and Keiding, N.（1993）. Statistical Methods Based on Counting Processes, Springer-Verlag, New York.
Andersen, P. K., Esbjerg, S. and Serensen, T. I. A.（2000）. Multi-state models for bleeding episodes
and mortality in liver cirrhosis, Statistics in Medicine, 19, 587–599.
. Analysis of survival data by the proportional odds model, Statistics in Medicine,
Bennett, S.（1982）
2, 273–277.
Breslow, N. E.（1972）. Disucussion of Professor Cox’s paper, Journal of the Royal Statistical Society,
Series B, 34, 216–217.
. Statistical Methods in Cancer Research, Vol. I, IARC Scientiﬁc
Breslow, N. E. and Day, N. E.（1980）
Publications, No. 32, IARC, Lyon.
Bilias, Y., Gu, M. and Ying, Z.（1997）. Towards a general asymptotic theory for Cox model withstaggered entry, The Annals of Statistics, 25, 662–682.
. Linear regression with censored data, Biometrika, 66, 429–436.
Buckley, J. and James, I.（1979）. Semiparametric analysis of transformation models with cenChen, K., Jin, Z. and Ying, Z.（2002）
sored data, Biometrika, 89, 659–668.
Chen, Y. Q. and Jewell, N. P.（2001）. On a general class of semiparametric hazards regression models,
Biometrika, 88, 687–702.
. Inference for a semiparametric model with panel data, Biometrika,
Cheng, S. C. and Wei, L. J.（2000）
87, 89–97.
. Analysis of transformation models with censored data,
Cheng, S. C., Wei, L. J. and Ying, Z.（1995）
Biometrika, 82, 835–845.
Cheng, S. C., Wei, L. J. and Ying, Z.（1997）. Predicting survival probabilities with semiparametrictransformation models, Journal of the American Statistical Association, 92, 227–235.
Cox, D. R.（1972）. Regression models and life tables, Journal of the Royal Statistical Society, Series
B, 34, 187–220.
. Partial likelihood, Biometrika, 62, 269–276.
Cox, D. R.（1975）
Diggle, P. J., Heagety, P. H., Liang, K.-Y. and Zeger, S. L.（2002）. Analysis of Longitudinal Data,
2nd ed., Oxford, New York.
. On the linear transformation model for censored data,Fine, J. P., Ying, Z. and Wei, L. J.（1998）
Biometrika, 85, 980–986.
. Counting Processes and Survival Analysis, John Wiley
Fleming, T. R. and Harrington, D.（1991）& Son, New York.
Fujii, T., Kawahara, A., Basaki, Y., Hattori, S., Nakashima, K., Nakano, K., Shirouzu, K., Kohno, K.,
Yanagawa, T., Yamana, H., Nishio, K., Ono, M., Kuwano, M. and Kage, M.（2008）. Expres-

統計数理

134

第 57 巻

第 1 号 2009

sion of HER2 and estrogen receptor alpha depends upon nuclear localization of Y-box binding
protein-1 in human breast cancers, Cancer Research, 68, 1504–1512.
. A generalized Wilcoxon test for comparing arbitrary single-censored samples,
Gehan, E. A.（1965）
Biometrika, 52, 203–223.
Gill, R. D.（1980）. Censoring and Stochastic Integrals, Mathematical Center Tract 124, MathematicalCenter, Amsterdam.
. A simple test of the proportional hazards assumption,
Gill, R. D. and Schumacher, M.（1987）
Biometrika, 74, 289–300.
. Diagnostic plots to reveal functional
Grambsch, P. M., Therneau, T. M. and Fleming, T. R.（1995）
form for covariates in multiplicative intensity models, Biometrics, 51, 1469–1482.
. Checking the linear transformation models for clustered failure time observations,
Hattori, S.（2008a）
Lifetime Data Analysis, 14, 253–266.
Hattori, S.（2008b）. Regression diagnostics of the semiparametric proportional rate model for irregularly-spaced repeated measurements, Communications in Statistics, Theory and Methods（inpress）.
Hougaard, P.（1999）. Multi-state models: A review, Lifetime Data Analysis, 5, 239–264.
. Iterative partial least squares with right-censored data
Huang, J. and Harrington, D. P.（2005）
analysis: A comparison to other dimension reduction techniques, Biometrics, 61, 17–24.
Jin, Z., Lin, D. Y., Ying, Z. and Wei, L. J.（2001）. A simple resampling method by perturbing theminimand, Biometrika, 88, 381–390.
Jin, Z., Lin, D. Y., Wei, L. J. and Ying, Z.（2003）. Rank-based inference for the accelerated failure
time model, Biometrika, 90, 341–353.
. Rank regression analysis of multivariate failure time data
Jin, Z., Lin, D. Y. and Ying, Z.（2005）
based on marginal linear models, Scandinavian Journal of Statistics, 33, 1–23.
. On least-squares regression with censored data, Biometrika,
Jin, Z., Lin, D. Y. and Ying, Z.（2006）
93, 147–161.
. Omnibus tests of the martingale assumption in the analysis
Jones, C. L. and Harrington, D. P.（2001）
of recurrent failure time data, Lifetime Data Analysis, 7, 157–171.
Kalbﬂeisch, J. D. and Prentice, R.（2002）. The Statistical Analysis of Failure Time Data, 2nd ed.,John Wiley & Sons, New York.
. Goodness-of-ﬁt test for the additive risk model with
Kim, J., Song, M. S. and Lee, S.（1998）(p > 2) -dimensional time-invariant covariates, Lifetime Data Analysis, 4, 405–416.
. Large sample theory of a modiﬁed Buckley-James estimator for
Lai, T. L. and Ying, Z.（1991）
regression analysis with censored data, The Annals of Statistics, 19, 1370–1402.
Lai, T. L. and Ying, Z.（1992）. Linear rank statistics in regression analysis with censored or truncateddata, Journal of Multivariate Analysis, 19, 13–45.
Li, L. and Gui, J.（2004）. Partial Cox regression analysis for high-dimensional microarray gene ex, i208–i215.
pression data, Bioinformatics, 20（Suppl. 1）
Li, L. and Li, H.（2004）. Dimension reduction methods for microarrays with application to censoredsurvival data, Bioinformatics, 20, 3406–3412.
. Goodness-of-ﬁt analysis for the Cox regression model based on a class of parameter
Lin, D. Y.（1991）
estimators, Journal of the American Statistical Association, 86, 725–728.
Lin, D. Y. and Geyer, J.（1992）. Computational methods for semiparametric linear regression withcensored data, Journal of Computational and Graphical Statistics, 1, 77–90.
. Model checking techniques for parametric regression with
Lin, D. Y. and Spiekerman, C. F.（1996）
censored data, Scandinavian Journal of Statistics, 23, 157–177.
. Goodness-of-ﬁt tests for the general Cox regression model, Statistica
Lin, D. Y. and Wei, L. J.（1991）

生存時間解析におけるセミパラメトリック推測とその周辺

135

Sinica, 1, 1–17.
. Semiparametric analysis of the additive hazard risk model, BiometLin, D. Y. and Ying, Z.（1994）
rika, 81, 61–71.
. Semiparametric analysis of general additive-multiplicative hazard
Lin, D. Y. and Ying, Z.（1995）
models for counting processes, The Annals of Statistics, 23, 1712–1734.
. Additive hazards regression models for survival data, Proceedings
Lin, D. Y. and Ying, Z.（1997）of the First Seattle Symposium in Biostatistics: Survival Analysis（eds. D. Y. Lin and T. R.
Fleming）, 185–198 Lecture Note in Statistics, No. 123, Springer-Verlag, New York.
. Semiparametric and nonparametric regression analysis of longitudinal
Lin, D. Y. and Ying, Z.（2001）
data（with discussion）, Journal of the American Statistical Association, 96, 103–126.
. Semiparametric regression analysis of longitudinal data with inforLin, D. Y. and Ying, Z.（2003）
mative drop-outs, Biostatistics, 4, 385–398.
Lin, D. Y., Wei, L. J. and Ying, Z.（1993）. Checking the Cox model with cumulative sums ofmartingale-based residuals, Biometrika, 80, 557–572.
. Conﬁdence bands for survival curves under theLin, D. Y., Fleming, T. R. and Wei, L. J.（1994）
proportional hazards model, Biometrika, 81, 73–81.
. Accelerated failure time models for counting processes,
Lin, D. Y., Wei, L. J. and Ying, Z.（1998）
Biometrika, 85, 605–618.
Lin, D. Y., Wei, L. J., Yang, I. and Ying, Z.（2000）. Semiparametric regression for the mean and ratefunctions of recurrent events, Journal of the Royal Statistical Society, Series B, 62, 711–730.
Lin, D. Y., Wei, L. J. and Ying, Z.（2002）. Model-checking techniques based on cumulative residuals,
Biometrics, 58, 1–12.
. Marginal regression of multivariate event times based on linear transformation models,Lu, W.（2005）
Lifetime Data Analysis, 11, 389–404.
. Additive risk models for survival data with highMa, S., Kosorok, M. R. and Fine, J.（2006）
dimensional covariates, Biometrics, 62, 202–210.
. A semiparametric additive regression model for longituMartinussen, T. and Scheike, T. H.（1999）
dinal data, Biometrika ，86, 691–702.
Martinussen, T. and Scheike, T. H.（2001）. Sampling adjusted analysis of dynamic additive regression
models for longitudinal data, Scandinavian Journal of Statistics, 28, 303–323.
. Eﬃcient estimation in additive hazards regression with
Martinussen, T. and Scheike, T. H.（2002a）current status data, Biometrika, 89, 649–658.
. A ﬂexible additive multiplicative hazards model, BiometMartinussen, T. and Scheike, T. H.（2002b）
rika, 89, 283–298.
Martinussen, T., Scheike, T. H. and Skovgard, I. M.（2002）. Eﬃcient estimation of ﬁxed and timevarying covariate eﬀects in multiplicative intensity models, Scandinavian Journal of Statistics,
29, 57–74.
. Maximum likelihood estimation in
Murphy, S. A., Rossini, A. J. and van der Vaart, A. W.（1997）
the proportional odds model, Journal of the American Statistical Association, 92, 968–976.
Nguyen, D. V. and Rocke, D. M.（2002）. Partial least squares proportional hazard regression for
application to DNA microarray survival data, Bioinformatics, 18, 1625–1632.
. Goodness-of-ﬁt methods for generalized linear mixed models, BioPan, Z. and Lin, D. Y.（2005）
metrics, 61, 1000–1009.
. A resampling method based on pivotal estimating
Parzen, M. I., Wei, L. J. and Ying, Z.（1994）
functions, Biometrika, 81, 341–350.
. Empirical Processes: Theory and Applications, Institute of Mathematical StatisPollard, D.（1990）
tics, Hayward, California.

統計数理

136

第 57 巻

第 1 号 2009

. Linear rank tests with right censored data, Biometrika, 65, 167–179.
Prentice, R. L.（1978）
Robins, J. M.（1993）. Estimation of the time-dependent accelerated failure time model in the presenceof confounding factors, Biometrika, 79, 321–334.
. Correcting for non-compliance in randomized trials using rank
Robins, J. M. and Tsiatis, A.（1991）
preserving structural failure time models, Communications in Statistics, Theory and Methods,
79, 2609–2631.
Scharfstein, D. O., Tsiatis, A. A. and Gilbert, P. B.（1998）. Semiparametric eﬃcient estimation in thegeneralized odds-rate class of regression models for right-censored time-to-event data, Lifetime
Data Analysis, 4, 355–391.
Scheike, T. H. and Martinussen, T.（2004）. On eﬃcient estimation and tests of time-varying eﬀects
in the proportional hazards model, Scandinavian Journal of Statistics, 31, 51–62.
. Direct modelling of regression eﬀects for transition probabilities
Scheike, T. H. and Zhang, M.（2007）
in multistate models, Scandinavian Journal of Statistics, 34, 17–32.
. Additive hazards Markov regression models illustrated with bone
Shu, Y. and Klein, J. P.（2005）
marrow transplant data, Biometrika, 92, 283–301.
Spiekerman, C. F. and Lin, D. Y.（1998）. Marginal regression models for multivariate failure timedata, Journal of the American Statistical Association, 93, 1164–1175.
. Nonparametric model checks in censored regression, Communications in Statistics,
Stute, W.（2000）
Theory and Methods, 29, 1611–1629.
Stute, W., Manteiga, G. and Quindimil, M. P.（1998）. Bootstrap approximations in model checks forregression, Journal of the American Statistical Association, 93, 141–149.
. A lack-of-ﬁt test for the mean function in a generalized linear model,
Su, J. Q. and Wei, L. J.（1991）
Journal of the American Statistical Association, 86, 420–426.
Sun, Y. and Wu, H.（2005）. Semiparametric time-varying coeﬃcients regression models for longitudinal data, Scandinavian Journal of Statistics, 32, 21–47.
. Martingale-based residuals for survival models, BioTherneau, T. M. and Grambsch, P. M.（1990）
metrika, 77, 147–160.
. A large sample study of Cox’s regression model, The Annals of Statistics, 9,Tsiatis, A. A.（1981）
93–108.
. Estimating regression parameters using linear rank tests for censored data,
Tsiatis, A. A.（1990）
The Annals of Statistics, 18, 354–372.
van der Vaart, A. W. and Wellner, J. A.（1996）. Weak Convergence and Empirical Processes,
Springer-Verlag, New York.
. Nonparametric pathway-based regression models for analysis of genomic
Wei, Z. and Li, H.（2007）
data, Biostatistics, 8, 265–284.
Wei, L. J., Ying, Z. and Lin, D. Y.（1990）. Linear regression analysis of censored survival data basedon rank tests, Biometrika, 77, 845–851.
. Maximum likelihood estimation of misspeciﬁed models, Econometrica, 50, 1–25.
White, H.（1982）
Yang, S. and Prentice, R.（1999）. Semiparametric inference in the proportinal odds regression model,
Journal of the American Statistical Association, 94, 125–136.
. Additive hazards model with multivariate failure time data, Biometrika,
Yin, G. and Cai, J.（2004）
91, 801–818.
Ying, Z.（1993）. A large sample study of rank estimation for censored regression data, The Annalsof Statistics, 21, 76–99.
. A test of ﬁt for a semiparametric additive risk model,
Yuen, K. C. and Burke, M. D.（1997）
Biometrika, 84, 631–639.
. Eﬃcient estimation for the accelerated failure time model, Journal
Zeng, D. and Lin, D. Y.（2007a）

生存時間解析におけるセミパラメトリック推測とその周辺

137of the American Statistical Association, 102, 1387–1395.
. Maximum likelihood estimation in semiparametric regression models
Zeng, D. and Lin, D. Y.（2007b）with censored data, Journal of the Royal Statistical Society, Series B, 69, 507–564.

138

Proceedings of the Institute of Statistical Mathematics Vol. 57, No. 1, 119–138(2009)

Semiparametric Inference in Survival Analysis and Related Topics
Satoshi Hattori
Biostatistics Center, Kurume University

Inference procedures for semiparametric non-proportional hazard regression models
such as an additive hazards model, a proportional odds model and an accelerated failure
time model have been developed. While the maximum partial likelihood method was
not applied successfully to these models, inference procedures have been developed forthem in a unified way based on martingale estimating equations. Furthermore, regression
diagnostics methods have been developed in a unified way based on martingale residuals.
This paper reviews recent developments of semiparametric failure time regression models
and some related topics.

Key words: Cox proportional hazards model, accelerated failure time model, additive hazards model,
empirical process, linear transformation model, martingale.