Pythonで生存時間解析 〜Nelson-Aalenによるハザード比の推定〜

12311 ワード
Python 生存時間解析 lifelines データ分析 統計学 Python テキストリンク

生存時間解析の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 として利用可能である。引き続き、政権データセットを用いる（※前記事参照）

In
T = data["duration"]
E = data["observed"]

from lifelines import NelsonAalenFitter
naf = NelsonAalenFitter()

naf.fit(T,event_observed=E)

モデルの適合後、NelsonAalenFitterは cumulative_hazard_ プロパティをDataFrameとして持っている。

In
print(naf.cumulative_hazard_.head())
naf.plot()

累積ハザード関数は生存関数ほど明確な考察を与えないが、ハザード関数は生存時間解析における先進的手法の基本である。累積ハザード関数$H(t)$を推定したことを思い出すと、この曲線の変化量がハザード関数の推定量であることが分かる。

In
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() 関数も用意されている。

In
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);

どちらのグループはより大きなハザードを持っているかより明確になり、非民主政権は一定のハザードを持っているように見える。

バンド幅を選択する明確な方法はなく、異なるバンド幅は異なる推論結果を導くため、選択には細心の注意が必要である。

In
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);

参考