【統計応用・医薬生物学】Cox比例ハザードモデルと尤度関数【統計検定1級対策】

生存時間解析の勉強を進めて、今回はCox比例ハザードモデルについて過去問に対応できるように知識をつけていきたいと思います。

ハザード関数と生存関数の知識が前提に必要なので、わからなかったらこちらをどうぞ。

medibook.hatenablog.com

目次：

Cox比例ハザードモデルとは？
生存時間解析における尤度関数を考える
Cox比例ハザードモデルにおけるβの推定

Cox比例ハザードモデルとは？

Cox比例ハザードモデルは回帰分析の一種で、生存率に関係してくる因子を説明変数、ハザード関数を目的変数として設定するモデルです。

例えば、あるi番目の被験者について、生存率に関係する3つの説明変数 $x_1i, x_2i, x_3i$ があるとします。それぞれに対応する偏回帰係数を $\beta_1, \beta_2, \beta_3$ とするとi番目の被験者のハザード関数は

$h_i(t)=exp(\beta_1x_1i+\beta_2x_2i+\beta_3x_3i)h_0(t)$

と表せます。

最後の $h_0(t)$ はベースラインハザード関数と呼ばれ、すべての説明変数が0のときのハザード関数を指します。

行列を使うとよりシンプルな形で式を構成できます（行列はなんとなく苦手ですが）。

上記の式のように、回帰係数の行ベクトルをβ、説明変数の列ベクトルをXiとおくと

$h_i(t)=exp(\beta X_i)h_0(t)$

と書けます。

この中の $exp(\beta X_i)$ がハザード比と呼ばれます。βが分かれば各説明変数がどれくらい生存率に影響を与えるかがわかるので、これが通常知りたい内容となります。

生存時間解析における尤度関数を考える

パラメータであるβを推定するには尤度関数の考え方が必要です。

2018年の統計応用の過去問でも出ています。

生存時間解析における尤度関数は各被験者の生存時間を $t_1, t_2,....,t_n$ とし、打ち切りがあった場合は（つまり観察期間中にイベントが発生しなかった場合） $\delta_i=0$ 、イベントが発生した場合は $\delta_i=1$ とするような定義関数を使って、以下の式で考えることができます。

$L(\beta)=\Pi_{i=1}^n\{f(t_i)\}^{\delta_i}\{S(t_i)\}^{1-\delta_i}$

イベントが発生した人にいては $f(t_i)$ の確率で発生しており、発生しなかった人は生存関数 $S(t_i)$ の確率で発生していないと言えるので、各被験者が独立であることを踏まえると、すべての積が尤度関数になると言えます。

Cox比例ハザードモデルにおけるβの推定

先程の尤度関数をCox比例ハザードモデルを使って書き換えてみます。

イベントが発生する確率 $f(t_i)=h_i(t)S(t_i)$ となり、 $S(t_i)=-logH(t_i)=exp(-\int_0^{t_i}h(t)dt)$ なので、まずはハザード関数を使って書き換えると

$L(\beta)=\Pi_{i=1}^n\{h_i(t)\}^{\delta_i}\{exp(-\int_0^{t_i}h_i(t)dt)\}$

となります。

次にβを推定するため、対数尤度関数を考えます。対数尤度＝0となるようなパラメータβを考えれば、βの推定を行うことができそうです。

対数尤度関数を考えると

$logL(\beta)=\sum\{\delta_ilogh(t_i)+log(exp(-\int_0^{t_i}h_i(t)dt))\}\\=\sum\{\delta_ilog(h_0(t_i))+\beta X_i\}-\sum\{exp\beta X_i\int_0^{t_i}h_0(t)dt\}$

となります。

ですが。

ここまでやっといて何なのですが、正攻法でβの推定をするならこのような計算になり、ベースラインハザード関数 $h_0(t)$ を求める必要が出てきます。

しかしながら実際問題興味があるのはそこではなくあくまでβの推定です。

そこでCox比例ハザードモデルでは上記の式における $h_0(t_i)$ の情報をまったく捨て去って、βを推定します。

あるj番目の被験者が、ある時点 $t_j$ でイベントが発生したときの確率は

$\frac{説明変数x_jをもつ被験者がイベント発生するハザード}{t_j時点で生存する被験者が1人イベント発生するハザード}$

と言えるので、式としては

$\frac{h_j(t_j)}{\sum h_l(t_j)}$

となります。分母はその時点で生存するすべての被験者（＝リスク集合）のハザードの総和となります。

これも先程のCox比例ハザードモデルの式を代入すると、ベースラインハザード関数が分母と分子で相殺されてβのみが残ります。

$\frac{h_j(t_j)}{\sum h_l(t_j)}=\frac{exp\beta X_j}{\sum exp\beta X_l}$

これをある時点 $t_j$ だけでなく、全ての時点について掛け合わせれば良いので、イベント発生した回数をrとすると、尤度関数は

$\Pi_{j=1}^r\frac{exp\beta X_j}{\sum exp\beta X_l}$

となります。ここからニュートンラフソン法という方法を用いて、βの最尤推定を行います。それ以上の詳細は出ないと思うので（あと自分も分からないので、、、）ここまでとします。

式を見てわかるように生存時間がどれだけであったか、というデータはバッサリと切り捨てられているため、この尤度は部分尤度関数と言われます。

この方法であればベースラインハザード関数を推定せずに済んでいるため、ベースラインハザード関数については何らかの確率分布を前提としていません。イベント発生する順番のみが推定に関係しており、イベント発生までの時間も結果に関係してこないため、この部分はノンパラメトリックな手法であると言えます。ただ、比例ハザードの部分は対数線形モデルとなっており、パラメトリックな手法であるため、二つを合わせてセミパラメトリックなモデルである、と言われます。

医学研究ではもっともよく用いられている手法なので、過去問でもほぼ毎回出題されているように見えますが、時間に関するデータがまったく入っていないというのは果たして本当に実際の状況を表せているのか不安になります。ある疾患の研究に対して何が本当に良い統計手法なのかっていう話は突き詰めると疑問だらけですね。

参考文献：

統計学入門−第11章

深堀りされすぎてもはや理解できない領域に達していますが、Cox比例ハザードモデルについて式含めて解説されています。

https://www.heisei-u.ac.jp/ba/fukui/pdf/analysis31.pdf

2018年の過去問に出ていた指数分布に従う場合の最尤推定が説明されています。ワイブル分布の最尤推定も載っています。

リンク

生存時間解析の本の中では式も丁寧でわかりやすいと思います。