コックス比例ハザードモデル

はじめに

コックス比例ハザード回帰は、生存率の推定値を調整して予測子変数の影響を定量化するためのセミパラメトリック法です。この手法は、説明変数の影響を一般的なベースラインハザード関数 h₀(t) の乗数として表します。ハザード関数はコックス比例ハザード回帰関数のノンパラメトリック部分ですが、予測子変数の影響は対数線形回帰です。0 に対するベースラインの場合、このモデルは次に対応します。

$h (X_{i}, t) = h_{0} (t) \exp [\sum_{j = 1}^{p} x_{i j} b_{j}],$

ここで、 $X_{i} = (x_{i 1}, x_{i 2}, \dots, x_{i p})$ は i 番目の被験者についての予測子変数、h(X_i,t) は時間 t における X_i についてのハザード率、h₀(t) はベースラインハザード率関数です。

コックス比例ハザードモデルをデータに当てはめるには、関数 coxphfit または関数 fitcox を使用します。fitcox の方が新しい関数であり、モデルに関する詳細な情報を含む CoxModel オブジェクトを返します。fitcox の使用例については、コックス比例ハザードモデルオブジェクトを参照してください。

ハザード率

コックス比例ハザードモデルは、値が X_i のときの個体またはアイテムのハザード率をベースライン値のときの個体またはアイテムのハザード率に関連付けます。このモデルは、次のようにハザード率の推定値を生成します。

$H R (X_{i}) = \frac{h (X_{i}, t)}{h_{0} (t)} = \exp [\sum_{j = 1}^{p} x_{i j} b_{j}] .$

このモデルは、ベースラインハザード関数が時間 t に依存しているが、予測子変数は時間に依存していないという仮定に基づいています。この仮定は比例ハザードの仮定とも呼ばれ、個体に対するハザード率が時間が経過しても変化しないことを示します。

ハザード率は、予測子変数値 X_i をもつ個体またはアイテムの瞬時故障の相対リスクを、ベースライン値をもつ個体またはアイテムと比較して表します。たとえば、予測子変数が喫煙ステータスであり、非喫煙がベースラインカテゴリである場合、ハザード率は、ベースラインカテゴリ (非喫煙者) と比較した喫煙者の相対的な瞬時故障率を示します。X^* に対するベースラインと予測子変数値 X_i の場合、ハザード率は次のようになります。

$H R (X_{i}) = \frac{h (X_{i}, t)}{h (X^{*}, t)} = \exp [\sum_{j = 1}^{p} (x_{i j} - x_{j}^{*}) b_{j}] .$

たとえば、ベースラインが予測子変数の平均値 (mean(X)) である場合、ハザード率は次のようになります。

$H R (X_{i}) = \frac{h (X_{i}, t)}{h (\bar{X}, t)} = \exp [\sum_{j = 1}^{p} (x_{i j} - {\bar{x}}_{j}) b_{j}] .$

ハザード率は生存率に関連しており、説明変数の値が X_i である個体の時刻 t における生存率は次のようになります。

$S_{X_{i}} (t) = S_{0} {(t)}^{H R (X_{i})},$

ここで S₀(t) は、ベースラインハザード率関数 h₀(t) を含む生存時間関数で、HR(X_i) は、ベースライン値に相対的な予測子変数値 X_i のハザード率です。

コックス比例ハザードモデルの拡張

比例ハザード (PH) の仮定を満たさない変数がある場合、層化コックスモデルおよび時間に依存する変数があるコックスモデルという 2 つのコックス比例ハザードモデルの拡張を使用することを検討できます。

PH の仮定を満たさない変数がカテゴリ化可能である場合は、層化コックスモデルを使用します。

$h_{s} (X_{i}, t) = h_{0 s} (t) \exp [\sum_{j = 1}^{p} x_{i j} b_{j}],$

ここで、添字 s は s 番目の階層を示します。層化コックスモデルでは、ベースラインハザード率関数が各階層で異なりますが、係数は共有されます。したがって、予測子変数の値が同じである場合、すべての階層でハザード率が同じになります。名前と値のペア 'Strata' を使用すると、階層化変数を coxphfit に含めることができます。コックスモデルオブジェクトでの層化コックスモデルの使用例については、コックス比例ハザードモデルオブジェクトを参照してください。

PH の仮定を満たさない変数が時間に依存する変数である場合は、時間に依存する変数があるコックスモデルを使用します。

$h (X_{i}, t) = h_{0} (t) \exp [\sum_{j = 1}^{p_{1}} x_{i j} b_{j} + \sum_{k = 1}^{p_{2}} x_{i k} (t) c_{k}],$

ここで、x_ij は時間に依存しない予測子の要素、x_ik(t) は時間に依存する予測子の要素です。時間に依存する変数を coxphfit に含める方法の例については、共変量が時間に依存するコックス比例ハザードモデルを参照してください。

部分尤度関数

各注釈変数の影響の推定値、つまり、各注釈変数の影響の推定ハザード率は、他のすべての変数が一定であることを前提にすると、exp(b) となります。ここで b はその変数の係数推定値です。係数推定値は、モデルの部分尤度関数を最大化することにより得られます。比例ハザード回帰モデルの部分尤度関数は、イベントの観測順序に基づきます。これは、各故障時間について推定された故障の部分尤度の積です。n 個の異なる故障時間に n 回の故障があり、 $t_{1} < t_{2} < \dots < t_{n}$ である場合、部分尤度は次のようになります。

$L = \frac{H R (X_{1})}{\sum_{j = 1}^{n} H R (X_{j})} \times \frac{H R (X_{2})}{\sum_{j = 2}^{n} H R (X_{j})} \times \cdot \cdot \cdot \times \frac{H R (X_{n})}{H R (X_{n})} = \prod_{i = 1}^{n} \frac{H R (X_{i})}{\sum_{j = i}^{n} H R (X_{j})} .$

リスク集合 R_i を使用すると、部分尤度を次のように書き直すことができます。

$L = \prod_{i = 1}^{n} \frac{H R (X_{i})}{\sum_{j \in R_{i}} H R (X_{j})},$

ここで R_i は、調査対象の被験者の中で i 番目の故障時間までイベントを経験しない被験者のインデックス集合です。

尤度比検定を使用すると、モデルへの項の追加の有意性を評価できます。1 つ目のモデルに p 予測子変数があり、2 つ目のモデルに p + r 予測子変数がある 2 つのモデルについて考えてみます。この場合、2 つのモデルを比較すると、–2*(L₁/L₂) には自由度 r (検定されている項の数) をもつカイ二乗分布が含まれます。

同時発生イベントの部分尤度関数

同時発生イベントがある場合、coxphfit は厳密な部分尤度を計算する代わりに、Breslow (既定) または Efron の方法を使用してモデルの部分尤度を近似します。厳密な部分尤度の計算では、同時発生イベント時間についてのリスク集合の順列全体が含まれるので、大量の計算が必要です。

最もシンプルな近似法は、Breslow の方法です。この方法では、各同時発生集合について同じ分母を使用します。

$L = \prod_{i = 1}^{d} \prod_{j \in D_{i}} \frac{H R (X_{j})}{\sum_{k \in R_{i}} H R (X_{k})},$

ここで、d は異なるイベント時間の個数、D_i はイベント時間が i 番目のイベント時間に等しいすべての被験者のインデックス集合です。

Efron の方法は Breslow の方法より正確ですが、シンプルです。この方法では、同時発生イベントの分母を次のように調整します。

$L = \prod_{i = 1}^{d} \prod_{j \in D_{i}} \frac{H R (X_{j})}{\sum_{k \in R_{i}} H R (X_{k}) - \frac{j - 1}{d_{i}} \sum_{k \in D_{i}} H R (X_{k})},$

ここで、d_i は D_i 内のインデックスの個数です。

例として、はじめの 2 つのイベントが同時発生、つまり t₁ = t₂ および $t_{2} < t_{3} < \dots < t_{n}$ であるとします。Breslow の方法では、はじめの 2 つの項の分母は同じです。

$L = \frac{H R (X_{1})}{\sum_{j = 1}^{n} H R (X_{j})} \times \frac{H R (X_{2})}{\sum_{j = 1}^{n} H R (X_{j})} \times \frac{H R (X_{3})}{\sum_{j = 3}^{n} H R (X_{j})} \times \frac{H R (X_{4})}{\sum_{j = 4}^{n} H R (X_{j})} \times \cdot \cdot \cdot \times \frac{H R (X_{n})}{H R (X_{n})} .$

Efron の方法では、2 番目の項の分母を調整します。

$L = \frac{H R (X_{1})}{\sum_{j = 1}^{n} H R (X_{j})} \times \frac{H R (X_{2})}{0.5 H R (X_{1}) + 0.5 H R (X_{2}) + \sum_{j = 3}^{n} H R (X_{j})} \times \frac{H R (X_{3})}{\sum_{j = 3}^{n} H R (X_{j})} \times \frac{H R (X_{4})}{\sum_{j = 4}^{n} H R (X_{j})} \times \cdot \cdot \cdot \times \frac{H R (X_{n}, t_{n})}{H R (X_{n}, t_{n})} .$

近似法は、coxphfit の名前と値のペア 'Ties' を使用して指定できます。

観測値の頻度または重み

コックス比例ハザードモデルには、観測値の頻度または重みを組み入れることができます。w_i が i 番目の観測値の重みであるとします。すると、重みがあるコックスモデルの部分尤度は次のようになります。

重みがある部分尤度

$L = \prod_{i = 1}^{n} \frac{H R_{w} (X_{i})}{\sum_{j \in R_{i}} w_{j} H R (X_{j})},$
ここで

$H R_{w} (X_{i}) = \exp [\sum_{j = 1}^{p} w_{j} x_{i j} b_{j}] .$
重みがある部分尤度と Breslow の方法

$L = \prod_{i = 1}^{d} \prod_{j \in D_{i}} \frac{H R_{w} (X_{j})}{{[\sum_{k \in R_{i}} w_{k} H R (X_{k})]}^{\frac{1}{d_{i}} \sum_{j \in D_{i}} w_{j}}}$
重みがある部分尤度と Efron の方法

$L = \prod_{i = 1}^{d} \prod_{j \in D_{i}} \frac{H R_{w} (X_{j})}{{[\sum_{k \in R_{i}} w_{k} H R (X_{k}) - \frac{j - 1}{d_{i}} \sum_{k \in D_{i}} w_{k} H R (X_{k})]}^{\frac{1}{d_{i}} \sum_{j \in D_{i}} w_{j}}}$

観測値の頻度または重みは、coxphfit の名前と値のペア 'Frequency' を使用して指定できます。

参照

[1] Cox, D. R., and D. Oakes. Analysis of Survival Data. London: Chapman & Hall, 1984.

[2] Lawless, J. F. Statistical Models and Methods for Lifetime Data. Hoboken, NJ: Wiley-Interscience, 2002.

[3] Kleinbaum, D. G., and M. Klein. Survival Analysis. Statistics for Biology and Health. 2nd edition. Springer, 2005.

[4] Klein, J. P., and M. L. Moeschberger. Survival Analysis. Statistics for Biology and Health. 2nd edition. Springer, 2003.

参考

ecdf | coxphfit | ksdensity

コックス比例ハザード モデル