生存時間分析について学習したことまとめ（指数分布からCox比例ハザードモデルの導入まで）

<script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/3.2.2/es5/tex-mml-chtml.min.js">
</script>
<script type="text/x-mathjax-config">
 MathJax.Hub.Config({
 tex2jax: {
 inlineMath: [['$', '$'] ],
 displayMath: [ ['$$','$$'], ["\\[","\\]"] ]
 }
 });
</script>
# この記事について
データサイエンス初学者の筆者が生存時間分析を学習するにあたり、学んだことを自分なりにかみ砕いて記載するものとなります。数学的に厳密な内容というより、お気持ちを伝えることに重きを置いた内容となります。

# 前準備
## 指数分布について
生存時間分析において、指数分布というものの理解が前提になります。初めて知った方は、以下が参考になるかと思います。

https://www.youtube.com/watch?v=4Y5otbAwGlc

## 指数分布の確率密度関数
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3890153/7e1bff2d-1202-2e17-4a54-9861152bf4e0.png)

**$x$を時間として、あるイベントが次に発生する確率密度が$y$で与えられます**。 確率密度関数なので、確率は区間の面積によって与えられます。$\lambda$は、単位時間あたりのイベント発生回数になります。以下ではこれを$f(x)$と表すことにします。

:::note info
確率密度関数がわからない、という方向けに、ざっくり解説します。
確率密度関数では、**面積＝確率**です。仮に$x$を時刻としたとき、時刻1と時刻2の間にイベントが発生する確率は、関数の下の部分の面積で与えられます。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3890153/437f6296-994a-8218-ae40-7735cf04e33e.jpeg">
:::

## 指数分布の累積密度関数
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3890153/aae833a3-dd39-74df-3f3d-d07b44e5ca1c.png)

先ほどの確率密度関数を積分したものになります。生存時間分析においてはこの関数が非常に重要で、**$x$を時間としたとき、その時点で対象がイベント発生済みである確率が$y$で与えられます**。 これは累積密度関数なので、確率は$y$の値によって求められます。先ほどと同様、$\lambda$は、単位時間あたりのイベント発生回数です。以下では、指数分布の累積密度関数を$F(x)$と表すことにします。

:::note info
累積密度関数がわからない、という方向けに、ざっくり解説します。
確率密度関数では、**$y$＝確率**です。仮に$x$を時刻としたとき、時刻3までにイベントが発生済みである確率が、$y$で与えられます。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3890153/d5e679da-270d-1b71-5bac-8a414f1f928d.jpeg">
:::

## お気持ち
議論を進める上で使うのは累積密度関数になるので、その意味だけ理解できれば一旦OKかと思います。具体例を用いて解説します。

例）機械Aは、10年間に平均1回故障するといます。さて、機械Aが3年以内に故障する確率はどのように求められるでしょうか？

この例では、
イベント ＝ 故障
$\lambda$ = 0.1（10年に1回 → 1年に0.1回）
となるので、累積密度関数は、以下で与えられます。

```math
F(x) = 1 - e^{-0.1x}
```
ここで、**$x$を時間としたとき、その時点で対象がイベント発生済みである確率が$y$で与えられます**ので、$x$に3を代入して確率を求めることができます。

```math
F(3) = 1 - e^{-0.1×3} ≒ 0.2592
```

A. 約25.9%

# 生存時間分析について
生存時間分析とは、**あるイベントが発生するまでの予想される期間を分析する統計学の一分野**です。イベントは何でもよいのですが、医療分野で生死をイベントとして考えられることが多いので、生存時間分析という名前がついているようです。

## 生存関数
「お気持ち」で示した例では機械Aが3年以内に故障する確率が約25.9%であると求めましたが、ここから「**約74.1%の確率では、3年後も故障していない**」ということが言えそうです。この約74.1%というのは、以下の式で与えられます。

```math
1 - F(3) = 1- ( 1 - e^{-0.1×3} ) ≒ 0.7408
```

ここでは機械Aが3年後故障していない確率を求めることができましたが、イベント＝死亡とすることで、生命が3年後死亡していない、すなわち生命が3年後生存している確率を求めることができます。このようなことから、$S(x) = 1 - F(x)$には**生存関数**という名前が付けられています。

```math
生存関数：
S(x) = 1 - F(x)
```
これを、$y = S(x)$としてグラフにすると、以下のような形になります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3890153/a357a3de-c3ff-015a-52f7-739a40bbc86b.png)

この曲線は、**カプランマイヤー生存曲線**と呼ばれます。

:::note warn
ここでは「打ち切り」という概念を無視しています。
:::

## Cox比例ハザードモデルの導入
生存関数では、$\lambda$が与えられているという前提のもと、時間に対する対象の生存率を求めました。しかし実際、$\lambda$は事前に与えられるものではありませんし、病状や性別、人種などの様々な要因によって決まってきます。ここでは、生存率に影響する共変量から適切な$\lambda$を推定するアプローチである、Cox比例ハザードモデルについて解説します。

### ハザード関数
まず、$\lambda$を一般化しましょう。「前置き」では単位時間あたりのイベント発生回数と説明しましたが、これは「**時刻$x$まで生存していたという条件のもと、次の瞬間に死亡する定量的な量**」と言い換えることができ、次のような条件付き確率密度関数で求めることができます。

```math
ハザード関数：h(x) = \frac{f(x)}{S(x)}
```

ここで$f(x)$が$\lambda$に従う指数分布の確率密度関数であるとき、$h(x) = \lambda$となります。また、導出は割愛しますが、ハザード関数を求めることができると、以下の形で生存関数を求めることができます。

```math
生存関数の一般式：S(x∣X)=exp(−\int_o^xh(t∣X)dt)
```

すなわち、ハザード関数を求めることができれば、生存関数を求めることができるのです。

## Cox比例ハザードモデル
比例ハザードモデルでは、対数ハザードが、共変量の線形和になっていることを仮定します。

```math
比例ハザードモデル：h(x) = e^{\alpha + \beta_1 x_1+ \cdots +\beta_p x_p}
```

これに、時間依存の初期条件を掛けたものが、Cox比例ハザードモデルになります。

```math
Cox比例ハザードモデル：h(x) = h_0(x)e^{\alpha + \beta_1 x_1+ \cdots +\beta_p x_p}
```

:::note info
なぜ対数ハザードで仮定するかという点について、疑問に持たれる方もいるかもしれません。理由をここで詳しく解説することは省略しますが、この仮定により比例ハザード性という、イベントが時間依存でないことを保証できるということが理由の一つにある、と筆者は認識しています。
:::

指数分布を考えると、時間に対する生存率を求めることができますが、実際に統計分析を行う際には様々な共変量を考える必要があり、Cox比例ハザードモデルはそのような場合に役に立ちます。モデルを与えるためには$\beta$を推定する必要があり、Cox比例ハザードモデルにおいては部分最尤という概念を用いて推定するようですが、この記事では記載を省略します。

以上