原理的にはハザードを目的変数にし、生存率に影響を与える因子を説明変数つまり共変数にして重回帰分析を行えば、それらの因子が生存率に与える影響を解析することができます。 しかしその場合はハザード関数λ(t)の具体的な姿を規定する必要があり、それは現実的には非常に困難です。 そこでロジスティック回帰分析と同じように、次のようなリンク関数を用いた重回帰型モデルを想定します。 (→10.1 ロジスティック回帰分析の原理 (3) 一般化線形モデル)
| η=ln{ | λ(t|x1,…,xp) ――――――― λ0(t) |
}=β0+β1x1+…+βjxj+…+βpxp+ε (j=1,…,p) |
| λ(t|x1,…,xp) ――――――― λ0(t) |
=exp(β0+β1x1+…+βjxj+…+βpxp+ε) | =exp(η)=HR |
基準ハザード関数λ0(t)は全ての共変数の値が0の時のハザード関数であり、ハザード関数λ(t|x1,…,xp)はどれか1つ以上の共変数の値が0ではない時のハザード関数です。 そしてこのモデルは、ハザード関数は共変数によっても時間によっても形を変えず、λ(t|x1,…,xp)とλ0(t)の比を対数変換した値つまり対数ハザード比と共変数の間に線形関係があると仮定したものです。 ハザード比を対数変換するのは、共変数との関係をより線形に近づけるためです。
このモデルは、共変数はハザード関数の形には影響を与えず値だけに影響を与える、つまり図11.4でいえばλ(t)のグラフを単に上下に平行移動させるだけであり、しかもその移動比率を対数変換した値は共変数の値に比例するという、現実にはほとんど有り得ないものです。 しかしこのモデルには、ハザード関数の具体的な姿を想定する必要がないという大きなメリットがあります。
そこで、一般的な多変量生命表解析ではこのモデルを想定します。 これを「比例ハザードモデル(proportional hazard model)」といい、このモデルに基づいた多変量生命表解析のことを、「Coxの比例ハザードモデルによる重回帰型生命表解析」といいます。 数学的には、この手法は第2節で説明したコックス・マンテルの検定を多変量に拡張した手法に相当します。
比例ハザードモデルでは、共変数と生存関数S(t|x1,…,xp)の間に次のような関係があります。
基準生存関数S0(t)は、原理的には全ての共変数の値が0の時の生存関数です。 しかしS0(t)を求めるのは難しいため、実際のデータでは共変数を無視した時の累積生存率曲線をカプラン・マイヤー法によって求め、それをS0(t)にします。 このS0(t)は全ての共変数が平均値の時の生存関数に相当するため、その時の対数ハザード比ηが0になるように定数β0を調整します。
上記の関係から、共変数に任意の値を入れた時の累積生存率曲線を理論的に計算し、グラフを描くことができます。 その理論的な累積生存率曲線は、図11.5の赤い累積生存率曲線と青い累積生存率曲線ように、t=0の時の値が1で、それ以後は基準生存関数S0(t)を平行移動したような形になります。
比例ハザードモデルは重回帰型ですから、偏回帰係数βjは、他の共変数が一定で共変数xjだけが「1」増加した時に、対数ハザード比がいくつ変化するかを表す値、つまり対数ハザード比の変化量になります。 対数ハザード比の変化量は、次のように共変数が「1」増加した時の対数ハザード比になりますから、それを指数変換して元のハザード比に戻すと、これは他の共変数が一定で共変数xjだけが「1」増加した時に、ハザードが何倍になるかを表す値になります。 このためこのハザード比は、他の共変数の影響を取り除いた補正ハザード比になります。
| ηj0=ln{ | λj0(t) ――― λ0(t) |
} |
| ηj1=ln{ | λj1(t) ――― λ0(t) |
} |
| ηj1-ηj0=ln{ | λj1(t) ――― λ0(t) |
}-ln{ | λj0(t) ――― λ0(t) |
}=ln{ | λj1(t) ――― λj0(t) |
}=βj |
| ∴ | λj1(t) ――― λj0(t) |
=exp(βj)=HRj |
比例ハザードモデルは一般化線形モデルの一種であり、対数ハザード比の回帰誤差が特殊な分布になります。 そこで回帰誤差が近似的に正規分布すると仮定して、重回帰分析と同じように最小2乗法を利用して回帰分析を行う方法も考えられますが、普通は最尤法を利用した繰り返し近似計算によって回帰分析を行います。 ただし最尤法を正確に適用するためには、基準ハザード関数を具体的に規定する必要があります。 そこで実際の計算では、基準ハザード関数を無視して最尤解を近似計算するかなりラフな手法を用います。
例えば表11.8のデータに比例ハザードモデルを当てはめ、最尤法を利用して解を求めると次のようになります。 (注1)
| y=ln{ | λ(t|x1,x2) ――――― λ0(t) |
}=b0+b1・x1+b2・x2=-0.420-0.669x1+0.735x2 |
| 症例番号 | 生存期間(転帰) | 生存数/観察数 | 累積生存率 | 累積生存率の標準誤差 |
|---|---|---|---|---|
| 1 | 1 | 69/70 | 0.986 | 0.014 |
| 2 | 2 | 68/69 | 0.971 | 0.020 |
| 3 | 2 | 67/68 | 0.957 | 0.024 |
| 4 | 3 | 66/67 | 0.943 | 0.028 |
| 5 | 3 | 65/66 | 0.929 | 0.031 |
| 6 | 3 | 64/65 | 0.914 | 0.033 |
| 37 | 3 | 63/64 | 0.9 | 0.036 |
| 7 | 4 | 62/63 | 0.886 | 0.038 |
| 8 | 4 | 61/62 | 0.871 | 0.040 |
| 9 | 4 | 60/61 | 0.857 | 0.042 |
| 38 | 4 | 59/60 | 0.843 | 0.043 |
| 10 | 5 | 58/59 | 0.829 | 0.045 |
| 11 | 5 | 57/58 | 0.814 | 0.046 |
| 12 | 5 | 56/57 | 0.8 | 0.048 |
| 13 | 5 | 55/56 | 0.786 | 0.049 |
| 39 | 5 | 54/55 | 0.771 | 0.050 |
| 40 | 5 | 53/54 | 0.757 | 0.051 |
| 14 | 6 | 52/53 | 0.743 | 0.052 |
| 41 | 7 | 51/52 | 0.729 | 0.053 |
| 15 | 8 | 50/51 | 0.714 | 0.054 |
| 16 | 8 | 49/50 | 0.7 | 0.055 |
| 17 | 9 | 48/49 | 0.686 | 0.055 |
| 42 | 9 | 47/48 | 0.671 | 0.056 |
| 43 | 10 | 46/47 | 0.657 | 0.057 |
| 44 | 10 | 45/46 | 0.643 | 0.057 |
| 45 | 11 | 44/45 | 0.629 | 0.058 |
| 18 | 12 | 43/44 | 0.614 | 0.058 |
| 19 | 12 | 42/43 | 0.6 | 0.059 |
| 20 | 12 | 41/42 | 0.586 | 0.059 |
| 21 | 12 | 40/41 | 0.571 | 0.059 |
| 22 | 13 | 39/40 | 0.557 | 0.059 |
| 46 | 13 | 38/39 | 0.543 | 0.060 |
| 47 | 14 | 37/38 | 0.529 | 0.060 |
| 23 | 16 | 36/37 | 0.514 | 0.060 |
| 48 | 18 | 35/36 | 0.5 | 0.060 |
| 49 | 18 | 34/35 | 0.486 | 0.060 |
| 50 | 19 | 33/34 | 0.471 | 0.060 |
| 51 | 19 | 32/33 | 0.457 | 0.060 |
| 52 | 21 | 31/32 | 0.443 | 0.059 |
| 53 | 23 | 30/31 | 0.429 | 0.059 |
| 54 | 25 | 29/30 | 0.414 | 0.059 |
| 55 | 26 + | (29/29) | 0.414 | 0.059 |
| 24 | 27 | 27/28 | 0.399 | 0.059 |
| 56 | 27 | 26/27 | 0.385 | 0.058 |
| 25 | 28 | 25/26 | 0.370 | 0.058 |
| 26 | 28 | 24/25 | 0.355 | 0.057 |
| 57 | 28 | 23/24 | 0.340 | 0.057 |
| 58 | 28 + | (23/23) | 0.340 | 0.057 |
| 59 | 30 | 21/22 | 0.325 | 0.056 |
| 27 | 31 | 20/21 | 0.309 | 0.056 |
| 60 | 32 | 19/20 | 0.294 | 0.055 |
| 28 | 32 + | (19/19) | 0.294 | 0.055 |
| 29 | 33 | 17/18 | 0.278 | 0.054 |
| 61 | 33 + | (17/17) | 0.278 | 0.054 |
| 30 | 34 | 15/16 | 0.260 | 0.054 |
| 31 | 35 + | (15/15) | 0.260 | 0.054 |
| 62 | 35 + | (14/14) | 0.260 | 0.054 |
| 32 | 36 + | (13/13) | 0.260 | 0.054 |
| 63 | 37 | 11/12 | 0.239 | 0.053 |
| 33 | 44 + | (11/11) | 0.239 | 0.053 |
| 64 | 49 | 9/10 | 0.215 | 0.053 |
| 65 | 52 + | (9/9) | 0.215 | 0.053 |
| 66 | 54 | 7/8 | 0.188 | 0.053 |
| 34 | 54 + | (7/7) | 0.188 | 0.053 |
| 35 | 55 | 5/6 | 0.157 | 0.052 |
| 67 | 56 | 4/5 | 0.125 | 0.050 |
| 36 | 56 + | (4/4) | 0.125 | 0.050 |
| 68 | 58 + | (3/3) | 0.125 | 0.050 |
| 69 | 59 + | (2/2) | 0.125 | 0.050 |
| 70 | 60 + | (1/1) | 0.125 | 0.050 |
この生命表中の生存期間と累積生存率をプロットしたものが累積生存率曲線になり、それが基準生存関数S0(t)になります。 そしてS0(t)を利用すれば、共変数が任意の値の時の理論的生存関数を求めることができます。 例えば重症度が軽症で、治療が無い時と有る時の理論的生存関数は次のようになり、それらをグラフ化すると図11.5のようになります。 なお図11.5のグラフでは、理論的生存関数には脱落例をプロットしてありません。
この場合のx1とx2のハザード比は、それぞれ次のようになります。 そしてあまり意味はありませんが、ロジティック回帰分析と同様に、最尤法による解が漸近的に正規分布するという性質を利用して、偏回帰係数が0かどうか、つまりハザード比が1かどうかの検定を行うことができます。
表11.8の重症度を無視し、治療の有無についてコックス・マンテルの検定を適用すると次のようになります。 (注2)
図11.6は、治療の有無別に実際の累積生存率曲線を実線で描き、そこに比例ハザードモデルを利用して求めた理論的生存関数を破線で重ねて描いたものです。 この図を見ると、重症度の影響を補正すると治療の有無の生存率の差が少し大きくなることがわかると思います。
比例ハザードモデルは現実にはほとんど有り得ないモデルであり、最尤解を求める計算も非常にラフな近似計算ですから、計算結果の信頼性はかなり低くなります。 したがって計算結果が医学的に納得できないものだったならば、共変数の組み合わせを変えてみたり、モデルを変えてみたりして、色々と検討する必要があります。 信頼性の低い計算結果に振り回されるのは賢明ではありません。
| 共変数 | 観測期間 | 転帰 | ||||
|---|---|---|---|---|---|---|
| x11 | … | x1j | … | x1p | t1 | d1 |
| : | : | : | : | : | ||
| xi1 | … | xij | … | xip | ti | di |
| : | : | : | : | : | ||
| xn1 | … | xnj | … | xnp | tn | dn |
※di:死亡例は1、脱落例は0のダミー変数
・比例ハザードモデル
回帰誤差εiが正規分布すると仮定せず、最尤法によってβの最尤推定値ｂを求めます。
この式を解くにはλ0(t)を規定する必要があります。 そこでコックスは、この式を直接用いず、λ0(t)を任意の負ではない局外関数としたままｂを近似的に推定する手法を開発しました。
まず、上記の尤度関数をλ0(t)を含む部分と含まない部分に分けます。 そしてλ0(t)を含む部分がｂの推定に与える影響は少ないとして、λ0(t)を含まない部分だけで尤度を計算します。 これを部分尤度(partial likelihood)といい、次のように表されます。
：死亡例だけを掛け合わせる
：リスク集合について合計する
実際のデータには同時死亡例または同時脱落例があるため、それを考慮したブレスロー・ペトの近似方法(Breslow-Peto approximation method)によって部分尤度を計算します。
：dt例のｘを合計したもの
この対数部分尤度関数にニュートン・ラプソン法を適用すると、次のようになります。 (→10.3 ロジスティック回帰分析の計算方法 (注2))
：リスク集合R(t)に関する、指数的な重みwiを付けたxjの平均値相当の値
：リスク集合R(t)に関する、指数的な重みwiを付けたxjとxlの共分散相当の値
またはCt=1(ブレスロー・ペトの近似)
偏回帰係数が0かどうかの検定、つまりハザード比が1かどうかの検定は、最尤推定値の漸近的正規性を利用したワルドの検定によって行います。
またロジスティック回帰分析と同様に、共変数がない時の尤度つまり偏回帰係数が全て0の時の尤度と、共変数がp個の時の尤度の比を利用した尤度比検定によって、共変数全体の回帰の検定を行うことができます。 (→10.3 ロジスティック回帰分析の計算方法 (2) 最尤法を利用する方法 (注2))
の時有意水準100α％で有意
偏回帰係数の初期値ｂ0は、全て0にするのが一般的な方法です。 しかし、次のような方法で求めることもできます。 前節で説明したように、ハザード関数λ(t)が時間とは無関係に一定だとすると、生存関数S(t)は指数関数になります。
この場合、個体の理論的な生存時間は∞です。 そこで現実的にはS(t)が非常に小さな値eになった時に死亡すると仮定すると、生存時間tとハザードモデルの間には次のような関係があります。
ln(t0)は基準ハザード関数の生存時間に相当するため、共変数を無視した時の、対数変換した生存時間の平均値になります。 このため上記のモデルは、対数変換した生存時間ln(t)を目的変数にした、定数項を少し修正した重回帰モデルになります。 したがって対数変換した生存時間を目的変数にした重回帰分析を行い、その時の偏回帰係数の符号を反対にしたものを初期値ｂ0にすることができます。
表11.8のデータについて、実際に計算してみましょう。 まず表11.8の観察期間を対数変換したものを目的変数にし、治療の有無と重症度を説明変数にした重回帰分析を行うと次のような結果になります。
この重回帰式における偏回帰係数の符号を反対にしたものを、比例ハザードモデルにおける偏回帰係数の初期値にします。 比例ハザードモデルの定数項は、全ての共変数に平均値を代入した時の対数ハザード比が0になるように調整します。 このため、偏回帰係数だけをニュートン・ラプソン法で推定します。
更新されたｂ1を用いて同様の計算を繰り返すと、3回目で値が収束します。 そしてこの偏回帰係数と、x1の平均値0.485714とx2の平均値1.01429から定数項b03を求めます。 またＨ3-1の対角要素を利用して、偏回帰係数の検定を行うことができます。
(dt1：t=tにおけるx1=1の死亡例数)
(nt1：t=tにおけるx1=1の総例数)
(U：コックス・マンテルの検定における分子)
にすると、
(nt0：t=tにおけるx1=0の総例数)
(I：コックス・マンテルの検定における分母)
(b：コックスのβ推定値)
(z：コックス・マンテルの検定統計量)
以上のように、この場合の回帰係数b1はコックス・マンテルの検定におけるコックスのβと一致し、回帰係数の検定は連続修正をしないコックス・マンテルの検定と一致します。 (→11.2 生存率の比較方法 (注1))
表11.8の治療の有無について、比例ハザードモデルを適用した結果と、コックス・マンテルの検定を適用した結果は次のようになります。 両者の結果はわずかに異なっていますが、これは両者の有限修正の有無と連続修正の有無が原因です。
第5節へ
webmaster@snap-tck.com
Copyleft (C) 2000-2010 SNAP(Sugimoto Norio Art Production)