(注1)　表11.3.1を一般化して、比例ハザードモデルを当てはめると次のようになります。

比例ハザードモデル：
η=β₀１_n + β₁ｘ₁ + … + β_jｘ_j + … + β_pｘ_p + ε=Ｘβ + ε=ｙ + ε 　　ｙ=Ｘβ
　　 　　 　　

回帰誤差ε_iが正規分布すると仮定せず、最尤法によってβの最尤推定値ｂを求めます。

死亡例の確率：f(t_i|ｘ_i)=λ(t_i|ｘ_i)S(t_i|ｘ_i) 　　脱落例の確率：S(t_i|ｘ_i)
全体の尤度：
対数尤度：

この式を解くにはλ₀(t)を規定する必要があります。 そこでコックスはこの式を直接用いず、λ₀(t)を任意の負ではない局外関数としたままｂを近似的に推定する手法を開発しました。 まず上記の尤度関数をλ₀(t)を含む部分と含まない部分に分けます。 そしてλ₀(t)を含む部分がｂの推定に与える影響は少ないと仮定して、λ₀(t)を含まない部分だけで尤度を計算します。 これを部分尤度(partial likelihood)といい、次のように表されます。

部分尤度：
：死亡例だけを掛け合わせる 　　：リスク集合について合計する
R(t)：リスク集合、t=tまで生存していた死亡例と脱落例の集合

実際のデータには同時死亡例または同時脱落例があるので、それを考慮したブレスロー・ペトの近似方法(Breslow-Peto approximation method)によって部分尤度を計算します。

部分尤度：
対数部分尤度：
d_t：t=tにおける同時死亡例数 　　：d_t例のｘを合計したもの 　　w_i=exp(β'ｘ_i)

この対数部分尤度関数にニュートン・ラプソン法を適用すると次のようになります。 (→10.3 ロジスティック回帰分析の計算方法 (注2))

　　

：リスク集合R(t)に関する、指数的な重みw_iを付けたx_jの平均値相当の値
Ｈ_k=∇ｇ_k

：リスク集合R(t)に関する、指数的な重みw_iを付けたx_jとx_lの共分散相当の値
または C_t=1(ブレスロー・ペトの近似)：有限修正 (→1.8 科学的研究の種類 (注1))
N_t：t=tにおける観察対象数(生存例数＋死亡例数) 　　s_t：t=tにおける生存例数
ｂ_k+1=ｂ_k-Ｈ_k^-1ｇ_k

偏回帰係数が0かどうかの検定つまりハザード比が1かどうかの検定は、最尤推定値の漸近的正規性を利用したワルドの検定によって行います。

V(ｂ)=-Ｈ_k^-1≒Ｉ_f^-1 　　E(-Ｈ_k)=Ｉ_f：情報行列 　　[Ｉ_f]_jj^-1：Ｉ_f^-1の第j対角要素
検定：＞χ²(1,α)の時、有意水準100α％で有意
推定：100(1-α)％信頼区間
→ 下限: 　上限：

またロジスティック回帰分析と同様に、共変数がない時の尤度つまり偏回帰係数が全て0の時の尤度と、共変数がp個の時の尤度の比を利用した尤度比検定によって共変数全体の回帰の検定を行うことができます。 (→10.3 ロジスティック回帰分析の計算方法 (注2))

対数尤度差：
尤度比検定：-2(L_BP(0) - L_BP(β))=χ_β²＞χ²(p,α)の時、有意水準100α％で有意

偏回帰係数の初期値ｂ₀は全て0にするのが一般的な方法です。 しかし次のような方法で求めることもできます。 第3節で説明したように、ハザード関数λ(t)が時間とは無関係に一定とすると生存関数S(t)は指数関数になります。 そしてその時、個体の理論的な生存時間は∞です。 そこで現実的にはS(t)が非常に小さな値eになった時に死亡すると仮定すると、生存時間tとハザードモデルの間には次のような関係があります。

λ(t)=λ(定数)　　S(t)=exp(-λt)=e　　ln(e)=-λt 　　

ln(t)=ln(t₀) - (β'ｘ_i + ε)

ln(t₀)は基準ハザード関数の生存時間に相当するため、共変数を無視した時の対数変換した生存時間の平均値になります。 そのため上記のモデルは対数変換した生存時間ln(t)を目的変数にし、切片を少し修正した重回帰モデルになります。 したがって対数変換した生存時間を目的変数にした重回帰分析を行い、その時の偏回帰係数の符号を反対にしたものを初期値ｂ₀にすることができます。

表11.3.1のデータについて実際に計算してみましょう。 まず観察期間を対数変換したものを目的変数にし、治療の有無と重症度を説明変数にした重回帰分析を行うと次のようになります。

この重回帰式における偏回帰係数の符号を反対にしたものを、比例ハザードモデルにおける偏回帰係数の初期値にします。 比例ハザードモデルの切片は、全ての共変数に平均値を代入した時の対数ハザード比が0になるように調整します。 そのため偏回帰係数だけをニュートン・ラプソン法で推測します。

初期値：
　　 　　

更新されたｂ₁を用いて同様の計算を繰り返すと、3回目で値が収束します。 そしてこの偏回帰係数と、x₁の平均値0.485714とx₂の平均値1.01429から切片b₀₃を求めます。 またＨ₃^-1の対角要素を利用して偏回帰係数の検定を行うことができます。

　　 　　 　　
ｂ₄=ｂ₃ - Ｈ₃^-1ｇ₃≒ｂ₃
b₀₃=0.669106×0.485714 - 0.734747×1.01429=-0.420249
比例ハザードモデル：y=-0.420 - 0.669x₁ + 0.735x₂
○偏回帰係数の検定と推定

ハザード比：HR₁=exp(-0.669)=0.512
ハザード比の95％信頼区間　下限：HR_1L=exp(-1.217)=0.296 　上限：HR_1U=exp(-0.121)=0.886

ハザード比：HR₂=exp(b₂)=exp(0.735)=2.085
ハザード比の95％信頼区間　下限：HR_2L=exp(0.372)=1.451 　上限：HR_2U=exp(1.097)=2.996
○全回帰の尤度比検定
このモデルの尤度：L_BP(β₁,β₂)=-191.79 　　偏回帰係数が全て0の時の尤度：L_BP(0,0)=-201.434
χ_β²=-2×(-201.434 + 191.79)=19.287(p=6.4838×10^-5)＞χ²(2,0.05)=5.991

(注2)　比例ハザードモデルにおいて、説明変数が1つだけで、しかもそれが0または1の値を取るダミー変数の場合について考えてみましょう。

回帰係数の初期値b₁₀=0とすると
w_i=exp(b₁₀x_i1)=exp(0・x_i1)=1 　　
(d_t1：t=tにおけるx₁=1の死亡例数)
(n_t1：t=tにおけるx₁=1の総例数)
(U：コックス・マンテルの検定における分子)

有限修正をにすると
(n_t0：t=tにおけるx₁=0の総例数)
(I：コックス・マンテルの検定における分母)
(b：コックスのβ推定値)
(χ_o²：コックス・マンテルの検定統計量)

以上のように、この場合の回帰係数b₁はコックス・マンテルの検定におけるコックスのβと一致し、回帰係数の検定は連続修正をしないコックス・マンテルの検定と一致します。 (→11.2 生存率の比較方法 (注1))

表11.3.1の治療の有無に比例ハザードモデルを適用し、b₁の初期値を0としてニュートン・ラプソン法を1回だけ行った結果と、連続修正をしないコックス・マンテルの検定を適用した結果は次のように一致します。 ただしニュートン・ラプソン法を収束するまで行うと、比例ハザードモデルの結果はわずかに異なったものになります。 部分尤度法による最尤法は苦し紛れの近似計算のため、最尤解から少しずれた解に収束してしまうのです。

(注3)　一般に回帰分析の予測精度を表す指標は寄与率です。 しかしロジスティック回帰分析や重回帰型生命表解析のように、最尤法を利用した回帰分析はたいてい寄与率を求めることができません。 そこで対数尤度を利用した擬似寄与率を求めて寄与率の代用にします。 ところがコックスの比例ハザードモデルによる重回帰型生命表解析は、部分尤度という苦し紛れの手を用いるので擬似寄与率を求めることができません。

そこで予後予測の精度を表す指標としてC統計量(C-index、Concordance index)という値が提唱されています(Harrell et al、1996)。 これはモデルから予測される生存時間と、実際の生存時間の大小関係がどの程度一致しているかを表すノンパラメトリックな指標であり、両者の順位一致係数に相当します。 C統計量は次のように定義されています。

t_i、t_j：実際の生存時間 (i,j=1,…,n、i≠j) 　　f(t_i)、f(t_j)：モデルから予測される生存時間
t_iとt_jの大小関係の判定方法は一般化ウィルコクソンの2標本検定と同じ (→11.2 生存率の比較方法 (注4))
・u_ij=1：{f(t_i)＞f(t_j) かつ t_i＞t_j} または {f(t_i)＜f(t_j) かつ t_i＜t_j}
　モデルから予測される生存時間の大小関係と実際の生存時間の大小関係が一致している
・u_ij=0.5：{f(t_i)=f(t_j) かつ t_i＞t_j} または {f(t_i)=f(t_j) かつ t_i＜t_j}
　モデルから予測される生存時間は同じだが実際の生存時間は異なっている
・u_ij=0：{f(t_i)＞f(t_j) かつ t_i＜t_j} または {f(t_i)＜f(t_j) かつ t_i＞t_j}
　モデルから予測される生存時間の大小関係と実際の生存時間の大小関係が反対
・計算から除外：t_i=t_jまたは判定不能
(n_u：判定できた比較回数)

上記のようにi番目のデータとj番目のデータを比較して一致スコアを付け、それを全てのデータについて合計した値がC統計量です。 コックスの比例ハザードモデルの場合、生存時間を予測する関数f(t_i)の代わりに対数ハザード比関数η(t_i)を用いると便利です。 ただし対数ハザード関数を用いると、対数ハザードの大小関係が生存時間の大小関係とは反対になるので注意が必要です。

C統計量をロジスティック回帰分析に適用することもできます。 その場合、生存時間の代わりにイベントの有無を一致の指標にします。 つまりモデルから予測されるイベント発生確率の大小関係と、実際のイベント発生率の大小関係(発生=1＞非発生=0)がどの程度一致しているかを表す値がC統計量になります。 すると上記の定義から、実際のイベント発生率が同じデータは計算から除外され、イベント発生例とイベント非発生例の間の比較結果だけを合計することになります。

そのためC統計量はイベント発生例とイベント非発生例の間で、モデルから予測されるイベント発生確率の大小関係を総当りで比較し、イベント発生例のイベント発生確率が大きい時に「一致(1)」として一致率を求めた値になります。 これはマン・ホイットニィのU検定におけるU値、つまり2群の間で値の大小関係を比較し、値が大きい方を勝ち(Upper)とした時の勝ち数を比較回数で割って勝率にしたものと同じ値になります。

そしてこの勝率は、ROC分析におけるROC曲線のAUC(曲線下面積)に相当します。 したがってロジスティック回帰分析にC統計量を適用すると、C統計量はイベント発生群とイベント非発生群のロジットスコア(対数オッズ比)についてROC曲線を描いた時のAUCに相当することになります。 (→9.2 群の判別と診断率 (注4))

寄与率はモデルから予測される値と実際の値の計量的な一致度を表します。 それに対してC統計量はモデルから予測される生存時間と実際の生存時間の順序が一致している程度を表します。 そのため寄与率が1(100％)なら予測値と実際の値がぴったり一致しているのに対して、C統計量が1になっても予測値と実際の値の順序が一致しているだけで、値までぴったり一致しているとは限りません。 そのため予測精度を表す指標としては寄与率ほど良い指標ではありません。 でもコックスの比例ハザードモデルは擬似寄与率を求めることができないので、致し方なく予測精度を検討するための参考として用いることがあるようです。

しかし第6節で説明するパラメトリック生命表解析は厳密な最尤法を用いるので、疑似寄与率を求めることができます。 そのため苦し紛れの部分最尤法に基づいた比例ハザードモデルを用い、しかもあまり良い指標ではないC統計量を用いるよりも、厳密な最尤法に基づいたパラメトリック生命表解析と疑似寄与率を用いる方が合理的でしょう。 (→11.6 パラメトリック生命表解析)

(注1)で求めた比例ハザードモデルについてC統計量を求めてみましょう。

このモデルから求めた対数ハザード比yをハザードスコアとすると、この値は生存時間と反比例する値になります。 そこでこの値を表11.3.1の全ての症例について計算し、C統計量を求めると次のようになります。

ちなみに、表11.3.1に指数分布モデルによるパラメトリック生命表解析を適用した時の疑似寄与率は次のようになります。 (→11.6 パラメトリック生命表解析 (注2))

多変量予測モデルの予測能力の増加を評価する指標としてNRI(Net Reclassification Improvement)とIDI(Integrated Discrimination Improvement)という値が提唱されていて、これらの指標を生存時間解析やロジスティック回帰分析で用いる時があります。 しかしこれらは判別分析のような後ろ向き研究用の指標のため、前向き研究用の手法である生存時間解析やロジスティック回帰分析で用いるのは非合理です。 さらにこれらの指標はオーソドックスな指標――例えば寄与率や偏回帰係数――と比べるとあまり良い指標ではありませんが、一応、紹介しておきます。 (→9.5 変数の選択 (注2))