前口上 | 目次 | 第1章 | 第2章 | 第3章 | 第4章 | 第5章 | 第6章 | 第7章 | 第8章 | 第9章 | 第10章 |
第11章 | 第12章 | 第13章 | 第14章 | 第15章 | 第16章 | 第17章 | 第18章 | 第19章 | 第20章 | 付録 |
1 | 2 | 3 | 4 | 5 | 6 |
第3節で説明したように生存率を左右するのはハザードであることから、ハザードを目的変数にし、それに影響を与える因子を説明変数にして重回帰分析を行えば、原理的にはそれらの因子が生存率に与える影響を解析することができます。 しかしその場合はハザード関数λ(t)の具体的な姿を規定する必要があり、それは現実には困難な場合が多いと思います。 そこでロジスティック回帰分析と同じように、次のようなリンク関数を用いた重回帰型モデルを想定します。 (→10.1 ロジスティック回帰分析の原理 (3)一般化線形モデル)
基準ハザード関数(baseline hazard function)λ0(t)は全ての共変数の値が0の時のハザード関数であり、ハザード関数λ(t|x1,…,xp)はどれか1つ以上の共変数の値が0ではない時のハザード関数です。 そしてこのモデルは、ハザード関数は共変数によっても時間によっても形を変えず、λ(t|x1,…,xp)とλ0(t)の比を対数変換した値つまり対数ハザード比と共変数の間に線形関係があると仮定したものです。 ハザード比を対数変換するのは共変数との関係をより線形に近づけるためです。
このモデルは「共変数はハザード関数の形には影響を与えず、値だけに影響を与える」という仮定に基づいたものです。 図11.3.1でいえば共変数はλ(t)のグラフを単に上下に平行移動させるだけであり、しかもその移動比を対数変換した値は共変数の値に比例するという、現実にはほとんど有り得ないモデルです。 しかしこのモデルにはハザード関数の具体的な姿を規定する必要がないという大きなメリットがあります。
これを比例ハザードモデル(proportional hazard model)といい、このモデルに基づいた多変量生命表解析のことをコックス(Cox)の比例ハザードモデルによる重回帰型生命表解析といいます。 数学的には、この手法は第2節で説明したコックス・マンテルの検定を多変量に拡張した手法に相当します。
比例ハザードモデルでは、共変数と生存関数S(t|x1,…,xp)の間に次のような関係があります。
基準生存関数(baseline survival function)S0(t)は、原理的には全ての共変数の値が0の時の生存関数です。 しかしS0(t)を求めるのは難しいので、実際のデータでは共変数を無視した時の累積生存率曲線をカプラン・マイヤー法によって求め、それをS0(t)にします。 このS0(t)は全ての共変数が平均値の時の生存関数に相当するので、その時の対数ハザード比ηが0になるように切片β0を調整します。
上記の関係から、共変数に任意の値を入れた時の累積生存率曲線を理論的に計算してグラフを描くことができます。 その理論的な累積生存率曲線は、図11.4.1の赤色の累積生存率曲線と青色の累積生存率曲線のようになります。 これらのグラフはt = 0の時の値が1で、それ以後は基準生存関数S0(t)を平行移動したような形になることがわかると思います。
比例ハザードモデルは重回帰型ですから、偏回帰係数βjは、他の共変数が一定で共変数xjだけが「1」増加した時に対数ハザード比がいくつ変化するかを表す値つまり対数ハザード比の変化量になります。 対数ハザード比の変化量は共変数が「1」増加した時の対数ハザード比になるので、それを指数変換して元のハザード比に戻すと、他の共変数が一定で共変数xjだけが「1」増加した時にハザードが相対的に何倍になるかを表す値になります。 そのためこのハザード比は他の共変数の影響を取り除いた補正ハザード比(調整ハザード比)になります。
比例ハザードモデルは一般化線形モデルの一種であり、対数ハザード比の回帰誤差が特殊な分布になります。 そこで回帰誤差が近似的に正規分布すると仮定して、重回帰分析と同じように最小2乗法を利用して回帰分析を行う方法が考えられます。 しかし通常はロジスティック回帰分析と同様に、最尤法を利用した繰り返し近似計算によって回帰分析を行います。 ただし最尤法を正確に適用するためには基準ハザード関数を具体的に規定する必要があります。 そこで実際の計算では基準ハザード関数を無視して最尤解を近似計算するかなりラフな手法を用います。
例えば表11.3.1のデータに比例ハザードモデルを当てはめ、最尤法を利用して解を求めると次のようになります。 (注1)
また観察期間と転帰のデータを用いて、カプラン・マイヤー法によって生命表を求めると次のようになります。
症例番号 | 生存期間(転帰) | 生存数/観察数 | 累積生存率 | 累積生存率の標準誤差 |
---|---|---|---|---|
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.4.1のようになります。 なお図11.4.1のグラフでは理論的生存関数には脱落例をプロットしてありません。
この場合のx1とx2のハザード比はそれぞれ次のようになります。 そしてロジティック回帰分析と同様に、最尤法による解が漸近的に正規分布するという性質を利用して偏回帰係数が0かどうかの検定、つまりハザード比が1かどうかの検定と推定を行うことができます。 ただし多変量生命表解析も記述統計学的手法なので推測統計学的手法である検定とは相性が悪く、ほとんどの場合は統計的仮説検定ではなく単なる有意性検定になります。
表11.3.1の重症度を無視し、治療の有無についてコックス・マンテルの検定を適用すると次のようになります。 (注2)
この結果と比例ハザードモデルによる結果を比べると、重症度の影響を補正すると治療の有無のハザード比が1からより離れる、つまり治療有と治療無のハザードの差がより大きくなり、生存率の差がより大きくなることがわかります。 図11.4.2は治療の有無別に実際の累積生存率曲線を太い実線で描き、そこに比例ハザードモデルを利用して求めた理論的生存関数を細い実線で重ねて描いたものです。 この図を見ると、重症度の影響を補正すると治療の有無の生存率の差が少し大きくなることがわかると思います。
比例ハザードモデルは現実にはほとんど有り得ないモデルであり、最尤解を求める計算も非常にラフな近似計算です。 しかもこの手法もノンパラメトリック手法なので死亡例の順番を利用するだけで死亡時間という重要な情報は利用しません。 そのため死亡例の時間間隔が異なっていても結果は変わらず、計算結果の信頼性はかなり低くなります。 したがって計算結果が医学的に納得できないものだったら、共変数の組み合わせを変えてみたり、パラメトリックモデルを利用してみたりして、色々と検討する必要があります。 信頼性の低い計算結果に振り回されるのは賢明ではありません。 (注3) (→11.6 パラメトリック生命表解析)
共変数 | 観測期間 | 転帰 | ||||
---|---|---|---|---|---|---|
x11 | … | x1j | … | x1p | t1 | d1 |
: | : | : | : | : | ||
xi1 | … | xij | … | xip | ti | di |
: | : | : | : | : | ||
xn1 | … | xnj | … | xnp | tn | dn |
回帰誤差εiが正規分布すると仮定せず、最尤法によってβの最尤推定値bを求めます。
この式を解くにはλ0(t)を規定する必要があります。 そこでコックスはこの式を直接用いず、λ0(t)を任意の負ではない局外関数としたままbを近似的に推定する手法を開発しました。 まず上記の尤度関数をλ0(t)を含む部分と含まない部分に分けます。 そしてλ0(t)を含む部分がbの推定に与える影響は少ないと仮定して、λ0(t)を含まない部分だけで尤度を計算します。 これを部分尤度(partial likelihood)といい、次のように表されます。
実際のデータには同時死亡例または同時脱落例があるので、それを考慮したブレスロー・ペトの近似方法(Breslow-Peto approximation method)によって部分尤度を計算します。
この対数部分尤度関数にニュートン・ラプソン法を適用すると次のようになります。 (→10.3 ロジスティック回帰分析の計算方法 (注2))
偏回帰係数が0かどうかの検定つまりハザード比が1かどうかの検定は、最尤推定値の漸近的正規性を利用したワルドの検定によって行います。
またロジスティック回帰分析と同様に、共変数がない時の尤度つまり偏回帰係数が全て0の時の尤度と、共変数がp個の時の尤度の比を利用した尤度比検定によって共変数全体の回帰の検定を行うことができます。 (→10.3 ロジスティック回帰分析の計算方法 (注2))
偏回帰係数の初期値b0は全て0にするのが一般的な方法です。 しかし次のような方法で求めることもできます。 第3節で説明したように、ハザード関数λ(t)が時間とは無関係に一定とすると生存関数S(t)は指数関数になります。 そしてその時、個体の理論的な生存時間は∞です。 そこで現実的にはS(t)が非常に小さな値eになった時に死亡すると仮定すると、生存時間tとハザードモデルの間には次のような関係があります。
ln(t0)は基準ハザード関数の生存時間に相当するので、共変数を無視した時の対数変換した生存時間の平均値になります。 そのため上記のモデルは対数変換した生存時間ln(t)を目的変数にし、切片を少し修正した重回帰モデルになります。 したがって対数変換した生存時間を目的変数にした重回帰分析を行い、その時の偏回帰係数の符号を反対にしたものを初期値b0にすることができます。
表11.3.1のデータについて実際に計算してみましょう。 まず観察期間を対数変換したものを目的変数にし、治療の有無と重症度を説明変数にした重回帰分析を行うと次のようになります。
この重回帰式における偏回帰係数の符号を反対にしたものを、比例ハザードモデルにおける偏回帰係数の初期値にします。 比例ハザードモデルの切片は、全ての共変数に平均値を代入した時の対数ハザード比が0になるように調整します。 そのため偏回帰係数だけをニュートン・ラプソン法で推測します。
更新されたb1を用いて同様の計算を繰り返すと、3回目で値が収束します。 そしてこの偏回帰係数と、x1の平均値0.485714とx2の平均値1.01429から切片b03を求めます。 またH3-1の対角要素を利用して偏回帰係数の検定を行うことができます。
以上のように、この場合の回帰係数b1はコックス・マンテルの検定におけるコックスのβと一致し、回帰係数の検定は連続修正をしないコックス・マンテルの検定と一致します。 (→11.2 生存率の比較方法 (注1))
表11.3.1の治療の有無に比例ハザードモデルを適用し、b1の初期値を0としてニュートン・ラプソン法を1回だけ行った結果と、連続修正をしないコックス・マンテルの検定を適用した結果は次のように一致します。 ただしニュートン・ラプソン法を収束するまで行うと、比例ハザードモデルの結果はわずかに異なったものになります。 部分尤度法による最尤法は苦し紛れの近似計算なので最尤解から少しずれた解に収束してしまうのです。
そこで予後予測の精度を表す指標としてC統計量(C-index、Concordance index)という値が提唱されています(Harrell et al、1996)。 これはモデルから予測される生存時間と実際の生存時間の大小関係がどの程度一致しているかを表すノンパラメトリックな指標であり、両者の順位一致係数に相当します。 C統計量は次のように定義されています。
上記のようにi番目のデータとj番目のデータを比較して一致スコアを付け、それを全てのデータについて合計した値がC統計量です。 コックスの比例ハザードモデルの場合、生存時間を予測する関数f(ti)の代わりに対数ハザード比関数η(ti)を用いると便利です。 ただし対数ハザード関数を用いると対数ハザードの大小関係が生存時間の大小関係とは反対になるので注意が必要です。
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))