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

生存時間解析のPythonパッケージである `lifelines` を用いて、Nelson-Aalenによるハザード比の推定を行います。本記事は以下の記事の和訳になります。

* [Estimating univariate models - LIFELINES](https://lifelines.readthedocs.io/en/latest/Survival%20analysis%20with%20lifelines.html)
    - [Estimating the survival function using Kaplan-Meier](https://lifelines.readthedocs.io/en/latest/Survival%20analysis%20with%20lifelines.html#estimating-the-survival-function-using-kaplan-meier) 
    - [Estimating hazard rates using Nelson-Aalen](https://lifelines.readthedocs.io/en/latest/Survival%20analysis%20with%20lifelines.html#estimating-hazard-rates-using-nelson-aalen) 【本記事の対象】
    - [Estimating cumulative hazards using parametric models](https://lifelines.readthedocs.io/en/latest/Survival%20analysis%20with%20lifelines.html#estimating-cumulative-hazards-using-parametric-models)
    - [Other types of censoring](https://lifelines.readthedocs.io/en/latest/Survival%20analysis%20with%20lifelines.html#estimating-cumulative-hazards-using-parametric-models)

以下の記事の続きです。

* [Pythonで生存時間解析 〜Kaplan-Meierによる生存関数の推定〜 - Qiita](https://qiita.com/roki18d/items/3d725333ff1616a8977c)

# 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` として利用可能である。引き続き、政権データセットを用いる（※[前記事](https://qiita.com/roki18d/items/3d725333ff1616a8977c)参照）

```python:In
T = data["duration"]
E = data["observed"]

from lifelines import NelsonAalenFitter
naf = NelsonAalenFitter()

naf.fit(T,event_observed=E)
```

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

```python:In
print(naf.cumulative_hazard_.head())
naf.plot()
```

![スクリーンショット 2020-12-02 18.26.13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/7956cc76-4f72-be98-b6cf-e8830a8ed339.png)


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

```python: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");
```

![スクリーンショット 2020-12-02 18.27.11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/00a6b792-9431-f37f-efbf-36a1678f3fb6.png)

変化量に注目すると、民主政権・非民主政権ともに一定のハザードを持っているが、民主政権の方がより大きなハザード定数を持っているといえる。

## ハザード関数の平滑化

累積ハザード関数の解釈は難しいが、一方でほとんどの生存時間解析は累積ハザード関数で行われており、これが推奨される。

累積ハザード関数の代わりに、より解釈の容易なハザード関数を導出することができるが、注意点がある。導出にはカーネル平滑器が関連しており、平滑化の量を制御するバンド幅 (`bandwidth`) パラメータを指定する必要がある。この機能は `smoothed_hazard_()` および `smoothed_hazard_confidence_intervals_()` メソッドに含まれる。

推定値と信頼区間をプロットするための `plot_hazard()` 関数も用意されている。

```python: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);
```

![スクリーンショット 2020-12-02 18.28.17.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/8f869192-7f10-d420-47b1-69d2fa0ad832.png)


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

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

```python: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);
```

![スクリーンショット 2020-12-02 18.29.18.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/342646/c662cbc8-4706-5902-9842-6ed9075659c3.png)


# 参考

* [Estimating univariate models - LIFELINES](https://lifelines.readthedocs.io/en/latest/Survival%20analysis%20with%20lifelines.html)
* [生存時間解析について – ノンパラメトリック - AI LAB](https://ai-lab.lapras.com/ml/non-parametric/)