Cox比例ハザード回帰モデルについて数式ありで、できるだけわかりやすくまとめる

ランダム化比較試験で非常によく出てくるCox比例ハザード回帰モデルですが、その概念って結構理解しづらいように思います。数式だけの説明だと分かりにくいし、なしだと何となく理解できてるのかもやもやする感じになる（あと欠点やモデルの問題点が分からない）ので、ほどほどに両方用いつつ、説明してみます。

数式きっついかもしれませんが、あった方がより正確に概念を理解できるように思います。毎度のことながら統計学を大学でやったりしたわけではないので、間違いがありましたら質問・ツッコミ等お願いいたします。

生存時間解析とハザードについて
ハザード関数と累積ハザード関数、生存関数の関係性
- ＜補足＞ハザード関数と累積ハザード関数、生存関数の式変形
Cox比例ハザード回帰モデルとは何なのか
比例ハザード性って本当に成り立つのか
まとめ

生存時間解析とハザードについて

まずそもそもCox比例ハザード回帰モデルの”ハザード”ってなんやねん、というところから始めます。

”ハザード”というのは瞬間死亡率と書かれたりもしますが、瞬間死亡率とは「ある時点から次の瞬間の時点までにおける死亡率」という意味になります。

具体的な例で考えてみますと、日単位で観察している研究があるとして、30日目には10人の生存者がいたけれど、31日目までに1人死亡した、と言う場合、1/10=10（%）が30日時点での瞬間死亡率と言えます。*1

なので、このハザードが算出される研究というのは「時間＋イベントの発生」というデータをもつ研究に限られます。こうしたデータの解析方法は生存時間解析と呼ばれます。「時間」という数量データを扱うわけですが、通常の連続変量を比較するようなt検定やマン・ホイットニーのU検定、回帰分析として、重回帰分析などは使えません。

なぜかといえば、実際は全ての被験者に対してイベント発生までの時間が測定できるわけではない＝打ち切り(censoring)があるからです。すべての被験者が試験期間を完遂したとしても最後までイベント発生しない被験者が存在する以上は、打ち切りの一部に含まれるため、打ち切り０の試験というのはまずありません。

図：打ち切りの例とイメージ

こんな感じで、いろいろな理由でイベントが発生しないことがあるわけですね。

そうなると普通のデータ解析の方法ではできないですし、また基本的にはパラメトリックなデータ（何らかの確率分布に従うこと）にはならないので、ノンパラメトリックな方法が必要となります。

次に、このハザードの解析方法を考えるうえで前提となるハザード関数と累積ハザード関数、生存関数について説明します。

ハザード関数と累積ハザード関数、生存関数の関係性

Cox比例ハザード回帰モデルを考える上で外せないのがハザード関数の概念です。なので、まずはそこを押さえます。

先ほど出てきた”ハザード”は当然ながら時間ごとに変化していきます。つまり、時間を変数とした関数になるので、これをハザード関数 h(t) と呼びます。

さらに、このハザード関数を積分したものを累積ハザード関数 H(t) と定義します。

そして、ある時点までに生存している確率の関数が生存関数 S(t) です。これは累積ハザード関数と関連があり、H(t)=-log S(t)となります。なぜこうなるかは数式を変形していくとわかるのですが、あとで補足します。数式アレルギーの人は「とりあえず生存関数から累積ハザード関数が出せるんだな」と思ってもらえればよいと思います。

要約すると、それぞれはこういう関係になります。

つまり、どれかの関数が分かればそれぞれすべて計算できるわけです。この中で実際の研究から観察されるのは生存関数の推測値であり、生存関数を推測する方法のひとつがカプランマイヤー法です。経験分布関数（赤池情報量規準の記事で少し触れました）といって実際に研究で観察された数値を用いて関数を推測することで、生存関数を導くことができます。ちなみにこのカプランマイヤー法で推定した数値をグラフにしたものが、よく論文でみられるカプランマイヤー曲線ですね。

さて、これでハザード関数について説明は終わりにして、本題のCox比例ハザード回帰モデルの話に移ります。

以前書いたカプランマイヤー法の記事の話はこちら、多分ネットで検索したらもう少し分かりやすいのがあると思います（汗

実際の論文から統計を学んでみる③-ログランク検定は何をしているのか-カプランマイヤー曲線 –

経験分布関数の話はこちら

カルバック-ライブラー情報量〜赤池情報量規準(AIC)までの概略をわかりやすく③【統計検定1級対策】

＜補足＞ハザード関数と累積ハザード関数、生存関数の式変形

なぜ、H(t)=-logS(t)になるのか、という話です。積分と微分、確率密度関数と分布関数の話が出ますが、数式が苦手な人は読み飛ばしてください。

最初に、ハザードの定義を数式で表現します。

①被験者がある時点tまで生存している条件下で

②極小の時間t+Δまでの間にイベントが発生する

③単位時間当たりの確率とするためにΔで割る

ということが言えるので、分布関数F(t)をｔまでの時点でイベントが発生する確率と定義するとハザード関数は式としてはこうなります。*2

h(t)=\frac{1}{\Delta}P(t\lt T\leq t+\Delta|T\gt t)\\=\frac{1}{\Delta}\frac{P(t\lt T\leq t+\Delta,T\gt t)}{P(T\gt t)}\\=\frac{1}{\Delta}\frac{F(t+\Delta)-F(t)}{1-F(t)}\\=\frac{F'(t)}{1-F(t)}

縦棒｜は条件付きの意味を表します。

また分布関数は $P(T\leq t)=F(t), P(T\gt t)=1-F(t)$ という変形ができます。分布関数がどういうものかについては、こちらもご参考ください。

12-1. 累積分布関数とは | 統計学の時間 | 統計WEB

実はこの辺は前書いた記事とほぼ一緒です。

実際の医学論文から統計を学んでみるⅡ②-イベント数/人年データをNNTに直す方法- – 脳内ライブラリアン

両辺を積分すると

$(左辺)=H(t)$

$(右辺)=\int_0^t\frac{F'(u)}{1-F(u)}du\\=-log(1-F(t))$

さて、ここで分母の1-F(t)は（全確率）ー（それまでにイベントが発症している確率）なので、生存関数S(t)と等しいことが分かります。

よって

を導くことができました。

Cox比例ハザード回帰モデルとは何なのか

ようやく本題のCox比例ハザード回帰モデルです。Cox比例ハザードと略されたりしますが、まず要するに「回帰分析である」ということが大事です。

回帰分析とは説明変数と呼ばれる因子（原因になるもののデータ）を用いて結果変数（結果となる数値）を予想する式を立てる方法のことを言います。Cox比例ハザード回帰モデルは先ほど述べたある個人の”ハザード関数”を予想するためのモデルなんです。

例えば、i番目の被験者のハザード関数を $h_i(t)$ とすると、Cox比例ハザード回帰モデルはこんな感じの式になります。*1

h_i(t)=h_0(t)exp(\beta_1x_{1i}+\beta_2x_{2i}+\beta_3x_{3i}+...\beta_px_{pi})

なんだこの式、、、ってなるかもしれませんが、落ち着いて式を見てみます。

まず $h_0(t)$ はベースラインハザード関数と呼ばれ、全ての説明変数が0の時のハザード関数になります。一番基本的で何もない場合のハザード関数ということになりますね。

続いて後ろにくっついている

$exp(\beta_1x_{1i}+\beta_2x_{2i}+\beta_3x_{3i}+...\beta_px_{pi})$ ですが、xは説明変数で、各被験者のデータになります。βは説明変数ごとにかけられる係数で実際の計測データから推測されて求まる数値です。expはネイピア数で（）でくっついているのは乗数ですね。肩に乗っけるとみづらいのでこのように表記します。

まず、なんでexpなんか出てくるのかということですが、乗数を使うことで０になることを避けるためです。ここが０になったりすると、 $h_i(t)=0$ と定数になってあり得ないこと（ハザードが０の被験者＝イベントを起こす確率が常に０）になってしまうので、この方が都合が良いわけです。

また、たいていの場合、説明変数は式のように複数用意されます。例えば結果変数が脳梗塞発症のハザード関数なら、年齢、高血圧/脂質異常/糖尿病の有無、性別、脳梗塞の既往、などなどです。 $x_{1i}$ は、i番目の被験者の1番目の説明変数と言う意味で、 $x_{2i}$ は、i番目の被験者の2番目の説明変数と言う意味になります。ちなみに説明変数の設定の仕方によって求まるβも変わってくるので、何を説明変数に入れるかが非常に重要です。

さて、ランダム化比較試験でこのモデルを使う時、我々が基本的に一番興味があるのはハザード比（HR; Hazard Ratio）です。式を使って実際にハザード比がどうなるのかみてみます。

まず、1番目の説明変数が１のとき介入群、０のときコントロール群のように設定してみます。「他の説明変数の条件は同じ」と仮定します。介入した被験者のハザード関数を $h_i(t)$ , コントロール群のハザード関数を $h_c(t)$ とすると

$h_i(t)=h_0(t)exp(\beta_1+\beta_2x_{2i}+\beta_3x_{3i}+…\beta_px_{pi})$

$h_c(t)=h_0(t)exp(0+\beta_2x_{2i}+\beta_3x_{3i}+…\beta_px_{pi})$

となります。expの最初の項が違うだけですね。

では次に、ハザード比をみてみると

となりますね。こうして、Cox比例ハザード回帰モデルから導き出した式を使って、いつも論文で見ているハザード比が確認できるわけです。このβは前述したように、実験結果から得られる推測値なので、95%信頼区間があり、そこからハザード比の信頼区間も得られます。

ここで大事なのは $exp(\beta_1)$ が時間を表す変数であるtを含んでいないことです。つまり、Cox比例ハザード回帰モデルでは時間によらずハザード比が一定なのです。このことを比例ハザード性の仮定といいます。ここが大切かつ解釈に注意が必要な点になります。次の項目でまた説明します。

ちなみにCox比例ハザード回帰モデルはセミパラメトリックモデルだ、と言われたりしますが、個人的にはずっとこの意味が分からなかったんですけれども、もう一度式を見つめなおしてみます。

最初のベースラインハザード関数 $h_0(t)$ は特定の分布に従うものではなく、ノンパラメトリックです。

また後ろの

exp(\beta_1x_{1i}+\beta_2x_{2i}+\beta_3x_{3i}+...\beta_px_{pi})

はβに規定される特定の分布に従っており、パラメトリックです。なので2つ合わせていることからセミパラメトリックモデルと呼ばれます。

比例ハザード性って本当に成り立つのか

さて、頻繁に用いられるCox比例ハザード回帰ですが、ちょっと冷静に考えてみるとこの比例ハザード性ってあり得るんでしょうか。介入群とコントロール群で常に病気の発症の比率が一定ってなんかおかしくないですか？

そんな疑問に答えてくれるのがJAMAの統計コーナー”JAMA Guide to Statistics and Methods”でたまたま見つけた記事でした。*3

ハザード比が経過した時間によって変動してしまうような比例ハザード性に沿わない例を３つほど、分かりやすく挙げています。

1. No immediate Effect

即座に効果が出ない場合です。心血管イベントに対するスタチンvsプラセボの試験が例に挙げられています。確かにスタチン処方で血管に作用しイベント抑制するには明らかに時間を要することが分かります。降圧剤、ＤＭ治療による心血管イベントなんかも同様であることが予想されますね。

2. Immediate and Delayed Effects in Opposite Directions

これは初期と後での効果が反対になっている場合です。例として大腸がん発症に対する内視鏡スクリーニング検査のありvsなしの試験が挙げられています。試験初期には当然ながらスクリーニング検査群で発見が増えるので大腸がん発症（というか診断）のハザードは上がりますが、試験後期には、スクリーニング検査群で前がん病変が発見され除去されるため、大腸がんの発症率は抑えられます。

3. Variation in Disease Susceptibility

疾患の起こしやすさが個人によって異なるため変動を来す場合です。例としては閉経後の女性へのエストロゲン＋プロゲステロン投与vsプラセボで冠動脈疾患の発症の差をみています。この場合、試験初期にはホルモン投与群で冠動脈疾患のハザードが大幅に上がる(HR1.8)のですが、後期にはむしろホルモン投与群の方が低くなります(HR0.70)。結局結果としてはHR1.24に落ち着きます。これは、冠動脈疾患をもともと起こしやすい素因あったホルモン投与の被験者群がホルモンの投与という原因によって早々に病気を起こし、あとの群はそもそもの素因が少ないため、あまり起こさなくなったという解釈ができます。

こう見てみると、介入に即効性があって、しかも患者層の疾患の起こしやすさが均質で変わらない場合でなければ比例ハザードなんて成り立たないということになります。じゃあ、このモデルが意味ないのかと言うとそういうわけではなく、最後の例をみても分かるように、あくまで求められるのは研究の全経過を平均化したようなハザード比ですよ、ということを知っておくことが大事だと思います。

患者さんに応用するとしても、カプランマイヤー曲線を参考にしつつ、実際の介入による機序やそこで起きていることを解釈して（上の例で挙げたように）用いることが必要です。例えば、誰もやらないでしょうけど、余命いくばくかもない人にNo immediate effectsな治療をするのはナンセンスなわけで、もとの研究の試験期間をみながらそこは判断することが要るのかなと思います。

比例ハザード性があるかどうか自体を仮説検定する、と言う方法もできなくはないですが、先ほどの記事では「検出力が不十分なので不要である」と断言されています。むしろ上記のようなブレを確認する指標としては95%信頼区間のほうが妥当なのでしょう。