生存時間解析のPythonパッケージである lifelines を用いて、Nelson-Aalenによるハザード比の推定を行います。本記事は以下の記事の和訳になります。
以下の記事の続きです。
Nelson-Aalenによるハザード比の推定
生存関数は生存データセットを要約・可視化するのに優れた方法であるが、唯一の手段ではない。母集団のハザード関数$h(t)$に関心がある時、残念ながらKM推定量に変換することができない。しかし、次で表される累積ハザード関数に対するノンパラメトリックな推定量が存在する。
$$H(t)=\int_{0}^{t}\lambda(z)dz$$
累積ハザード関数に対する次の推定量はNelson-Aalen推定量と呼ばれる。
$$\hat{H}(t)=\sum_{t_i \leq t} \frac{d_i}{n_i}$$
ここで$d_i$は時点$t_i$における死亡数、$n_i$は死亡の可能性のあるものの数である。lifelines では、この推定量は NelsonAalenFitter として利用可能である。引き続き、政権データセットを用いる（※前記事参照）
T = data["duration"]
E = data["observed"]
from lifelines import NelsonAalenFitter
naf = NelsonAalenFitter()
naf.fit(T,event_observed=E)
モデルの適合後、NelsonAalenFitterは cumulative_hazard_ プロパティをDataFrameとして持っている。
print(naf.cumulative_hazard_.head())
naf.plot()
累積ハザード関数は生存関数ほど明確な考察を与えないが、ハザード関数は生存時間解析における先進的手法の基本である。累積ハザード関数$H(t)$を推定したことを思い出すと、この曲線の変化量がハザード関数の推定量であることが分かる。
naf.fit(T[dem], event_observed=E[dem], label="Democratic Regimes")
ax = naf.plot(loc=slice(0, 20))
naf.fit(T[~dem], event_observed=E[~dem], label="Non-democratic Regimes")
naf.plot(ax=ax, loc=slice(0, 20))
plt.title("Cumulative hazard function of different global regimes");
変化量に注目すると、民主政権・非民主政権ともに一定のハザードを持っているが、民主政権の方がより大きなハザード定数を持っているといえる。
ハザード関数の平滑化
累積ハザード関数の解釈は難しいが、一方でほとんどの生存時間解析は累積ハザード関数で行われており、これが推奨される。
累積ハザード関数の代わりに、より解釈の容易なハザード関数を導出することができるが、注意点がある。導出にはカーネル平滑器が関連しており、平滑化の量を制御するバンド幅 (bandwidth) パラメータを指定する必要がある。この機能は smoothed_hazard_() および smoothed_hazard_confidence_intervals_() メソッドに含まれる。
推定値と信頼区間をプロットするための plot_hazard() 関数も用意されている。
bandwidth = 3.
naf.fit(T[dem], event_observed=E[dem], label="Democratic Regimes")
ax = naf.plot_hazard(bandwidth=bandwidth)
naf.fit(T[~dem], event_observed=E[~dem], label="Non-democratic Regimes")
naf.plot_hazard(ax=ax, bandwidth=bandwidth)
plt.title("Hazard function of different global regimes | bandwidth=%.1f" % bandwidth);
plt.ylim(0, 0.4)
plt.xlim(0, 25);
どちらのグループはより大きなハザードを持っているかより明確になり、非民主政権は一定のハザードを持っているように見える。
バンド幅を選択する明確な方法はなく、異なるバンド幅は異なる推論結果を導くため、選択には細心の注意が必要である。
bandwidth = 8.0
naf.fit(T[dem], event_observed=E[dem], label="Democratic Regimes")
ax = naf.plot_hazard(bandwidth=bandwidth)
naf.fit(T[~dem], event_observed=E[~dem], label="Non-democratic Regimes")
naf.plot_hazard(ax=ax, bandwidth=bandwidth)
plt.title("Hazard function of different global regimes | bandwidth=%.1f" % bandwidth);
