生存分析

生存分析（survival analysis） は、あるイベントの発生確率や、イベント発生までの期間に関わる要因を調べたいときなどに使う。

「イベント」は分野や問題設定によって様々で、例えば医学・疫学においては「死亡」というイベントについて扱う。マーケティングにおいて顧客の「解約」をイベントとしたり、メーカーが機械の「故障」というイベントについて分析・予測するといった応用研究もある。

通常、イベントが発生した場合、それ以降その個体のデータは得られなくなる（脱落する；打ち切りデータ censored data となる）

（出所：生存時間分析入門 | 株式会社ARISE analytics（アライズアナリティクス））

データの例¶

import pandas as pd

df = pd.DataFrame([
    (True, 1),
    (True, 3),
    (False, 5),
    (True, 7),
    (True, 9),
    (False, 11),
    (False, 13),
    (True, 15),
], columns=["status", "days"])
df = df.sample(frac=1, random_state=0).reset_index(drop=True)
df["id"] = range(1, df.shape[0]+1)
df = df[["id", "status", "days"]]
df

Loading...

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
for status in df["status"].unique():
    idx = df["status"] == status
    color = "dimgray" if status else "white"
    ax.scatter(df.loc[idx, "days"], df.loc[idx, "id"], label=f"status={status}", color=color, edgecolors="dimgray")
    for _, row in df.loc[idx, :].iterrows():
        x1 = 0
        x2 = row["days"]
        y1 = y2 = row["id"]
        ax.plot([x1, x2], [y1, y2], color="dimgray")
ax.legend()
ax.set(xlabel="Survival Time (days)", ylabel="ID")
fig.show()

ハザード関数¶

生存率関数¶

ある患者の生存時間を表す確率変数を $T$ とする（ $T>0$ ）。

$T$ の確率密度関数（PDF）を $f(t)$ とする。 $f(t)$ は「期間tが経過した時点で死亡する確率」を意味する。

$T$ の累積分布関数（CDF）は

F(t) = P(T \leq t) = \int^t_0 f(u) du

(1)

$F(t)$ は「期間tが経過するまでに死亡している確率」である。

$F(t)$ をベースに、「期間tが経過するまで死亡していない確率」を意味する 生存率関数（survival rate function） あるいは 生存関数（survival function） と呼ばれる関数 $S(t)$ を考えることができる

S(t) = P(T > t) = 1 - F(t) = \int_t^\infty f(u) du

(2)

f(t) = - \frac{d S(t)}{d t}

(3)

という関係になる

生存率関数の推定¶

もし正確な生存時間が観測できているなら（研究終了時までに全固体が死亡しているなら）、生存関数は単純に

\hat{S}(t) = \frac{ \text{$t$日を超えて生存している個体数} }{ \text{全個体数} }

(4)

で推定できる。

例えば

	Status	Survival_in_days
1	True	8.0
2	True	10.0
3	True	20.0
4	False	25.0
5	True	59.0

というデータの場合、 $\hat{S}(t=11)$ は全員が11日目を観測しているので計算できるが、 $\hat{S}(t=30)$ は4番目の個体がまだ30日目を観測していないのでわからない

この場合、カプラン・マイヤー推定量（Kaplan–Meier estimator - Wikipedia）などで推定を行うことになる

生存曲線（カプラン・マイヤー曲線）¶

横軸に時間、縦軸に生存率をとった折れ線グラフ

例えば、次のようなデータがあったとする。

Loading...

import matplotlib.pyplot as plt
from sksurv.nonparametric import kaplan_meier_estimator

time, survival_prob, conf_int = kaplan_meier_estimator(
    event=data_y["Status"],
    time_exit=data_y["Survival_in_days"],
    conf_type="log-log"
)
plt.step(time, survival_prob, where="post")
plt.fill_between(time, conf_int[0], conf_int[1], alpha=0.25, step="post")
plt.ylim(0, 1)
plt.ylabel(r"est. probability of survival $\hat{S}(t)$")
plt.xlabel("time $t$")
plt.show()