下記の定義で与えられるハザード関数（単にハザードともいう）は、各時点における瞬間的な死亡のリスクや危険度を表すために使われます。
$$ h(t)=\lim_{\delta t \to 0}\frac{P(t\leq T<t+\delta t|T\geq t)}{\delta t} $$
生存関数同様に、ハザード関数に関しても生存時間データから近似的な値を推定することが可能です。 今回はそんなハザード関数の推定に関して見ていきます。
生命表法/保険数理法によるハザード関数の推定
生存関数を推定したときと同様に、生存時間データが $m$ 個のグループに分かれているとします。 単位時間当たりの平均的な死亡リスクは、その区間での死亡数を同じ区間内での平均生存時間で割ることで算出できます。
$j$ 番目の区間における死亡数を $d_j$、リスク集合の大きさを $n_j$、$j$番目の区間の長さを $\tau_j$ とし、区間内の死亡率が一定であると仮定すると、
$$ h^*(t) = \frac{d_j}{(n_j - d_j/2)\tau_j} $$ という関数で与えられます。
カプラン・マイヤー型のハザード関数の推定
区間でグループ化されていない生存時間データを用いる場合、カプラン・マイヤー型のハザード関数を用いることでハザード関数を推定することができます。
$j$ 番目の死亡時点 $t_{(j)}$ における死亡数を $d_j$、リスク集合の大きさを $n_j$ とし、$\tau_j = t_{(j+1)} - t_{(j)}$ とします。区間 $t_{(j)}$ から $t_{(j+1)}$ までの区間におけるハザード関数は、
$$ \hat{h}(t) = \frac{d_j}{n_j\tau_j} $$
によって推定することができます。
カーネル平滑化（カーネル密度推定）
引用：Kernel density estimation - Wikipedia
カーネル平滑化は、回帰関数や確率変数の確率密度関数の推定に使用するノンパラメトリック手法のひとつです。
ある母集団の標本のデータが与えられたとき、カーネル平滑化を行うことでその母集団のデータを外挿することができます。
$x_1, x_2, ..., x_n$ を（未知の）確率密度関数 $f$ を持つ独立同分布からの標本とします。カーネル関数 $K$、バンド幅（平滑化パラメータ）$b$ のとき、カーネル密度推定量は下記で表されます。
$$ \hat{f}(x, b) = \frac{1}{nb}\sum^{n}_{j=1}K\Big{(}\frac{x-x_j}{h}\Big{)} $$
ヒストグラムからカーネル密度推定を行う場合など、ここで用いられるカーネル関数は一般的に標準的なガウス関数（平均0、分散1）が用いられます。
参考：Kernel (statistics) - Wikipedia)
ハザード関数の平滑化
ハザード関数の推定値を時系列にプロットすると、凹凸が見られます。
この関数に上記で紹介した平滑化を施すことでパターンがよりはっきり見えるように関数を変換できます。
生存時間分析では、Epanechnikovのカーネル関数（$K(u) = 0.75(1-u^2) (-1 \leq u \leq 1)$）を用いたカーネル平滑化が行われます。
$$ h^\dagger(t) = \frac{1}{b}\sum^r_{j=1} 0.75 \bigg\{1 - \Big(\frac{t-t_{(j)}}{b}\Big)^2\bigg\} \frac{d_j}{n_j} $$
因みにこのカーネル関数はlilinesではutils/__init__.py内で下記のように定義されており、ハザード関数の推定値を算出するメソッド内（smoothed_hazard_）で呼ばれています。
def epanechnikov_kernel(t, T, bandwidth=1.0): M = 0.75 * (1 - ((t - T) / bandwidth) ** 2) M[abs((t - T)) >= bandwidth] = 0 return M
lifelinesでハザード関数の推定
lifelinesでは平滑化済みのハザード関数の値を算出する関数が実装されています。
from lifelines import NelsonAalenFitter naf = NelsonAalenFitter() naf.fit(durations=T, event_observed=E) naf.smoothed_hazard_(bandwidth=10)
この値を信頼区間とともにプロットする際はplot_hazard関数を使います。
naf.plot_hazard(bandwidth=10)
因みに、ハザード関数の生の値はクラス内に保持しない仕様になっています。 取得したい場合は、累積ハザード関数の差分をとる必要があります。
# ハザード関数の値を累積ハザード関数の差分から取得してプロット hazard_ = naf.cumulative_hazard_.diff().fillna(naf.cumulative_hazard_.iloc[0]) hazard_.plot(drawstyle="steps-post")
累積ハザード関数の推定
で見たように、時間 $t$ における累積ハザード関数 $H(t)$ は、生存時間を用いて表されることから、 生存時間データのモデリングには累積ハザード関数が重要となります。
$$ H(t) =-{\log{S(t)}} $$ 生存関数のカプラン・マイヤー推定量を $\hat{S(t)}$ とすると、累積ハザード関数は
$$ \hat{H}(t) = -{\log{\hat{S}(t)}} $$
によって推定することができます。
生存関数のカプラン・マイヤー推定量は下記であらわされるため、
$$ \hat{S}(t) = \prod_{j=1}^k\frac{n_j-d_j}{n_j} $$
累積ハザード関数の推定値は、
$$ \begin{align} \hat{H} (t) &= -{\log{\hat{S}(t)}} \\ &= -\sum_{j=1}^k{\log{(\frac{n_j-d_j}{n_j}})} \end{align} $$
となります。
なお、生存関数にネルソン・アーレン推定量を用いる場合、 生存関数のネルソン・アーレン推定量は下記であらわされるため、
$$ \hat{S}(t) = \prod_{j=1}^k\exp{(-d_j/n_j)} $$
累積ハザード関数の推定値は、
$$ \begin{align} \tilde{H}(t) &= -{\log{\tilde{S}(t)}} \\ &= \sum_{j=1}^k{d_j/n_j} \end{align} $$
となります。
lifelinesでは累積ハザード関数のネルソン・アーレン推定量を算出することができます。
累積ハザード関数のネルソン・アーレン推定量に関しては下記の記事で算出しているので参考にしてみてください。
www.kimoton.com
