【lifelines中文wiki】生存分析简介

【lifelines中文wiki】生存分析简介

英文原版连接

应用领域

生存分析最初是用来预测个人的寿命的。当遇到“这群人会活多长？”这类问题的时候，保险精算师、健康专家就会采用生存分析方法来解答。“这群人可以是一个国家的全体国民、或者患有某种疾病的人群、或者其它群体。一般来说，生存分析的研究对象是某一类患者。
生存分析的应用不局限于传统的“出生和死亡问题”，可以是用于预测任意事件的持续时间。例如，医学专业人士或许对两次“生孩子”之间的时间长度感兴趣，那么在这个问题中，“出生”指的是分娩，“死亡”指的是再次怀孕（可见，我们在放宽“出生和死亡”的定义）。再如，用户订购某种服务（例如订阅报纸），在这里，“出生”指的是“用户订购服务”这一事件，而“死亡”指的是“用户停止继续购买服务”这一事件。

数据删失

当你想要预测个体的生存时间时，往往会发现，并不是所有个体的“死亡”事件都已经发生。举个栗子，医学研究者不大可能等几十年，直到所有被观察者都死亡，再结束观察继而总结规律并做预测。他/她们对仅在几年或几个月后改善生命的有效性感兴趣。

我们习惯将生存分析的数据集中没有发生死亡事件的样本称为右删失样本——也就是说由于某些外部原因，直到数据收集结束那一刻，我们都没有观察到该样本死亡。对于右删失样本，我们唯一知道的是右删失样本至少活到了我们停止观察的那一刻。

注意：也有左删失样本——即我们没有观察到其“出生”时间。

数据分析中，我们经常犯的一个错误是：简单粗地选择忽略右删失样本。以下解释为啥不能忽略。假设有一群病患，该人群由A和B这两个子群组成。其中A这群人的生存时间很短，平均2个月。B这群病患的生存时间则长一些，平均12个月。但我们可能事先不知道这个区别。假设在第10个月，我们想要调查这群病患的平均生存时间。
下图中，红线表示在第10个月之前，死亡的病患，蓝线表示右删失样本（直到第10个月，仍然活着的病患）。要在第10个月估计这群病患的平均生存时间。如果我们简单粗暴地忽略右删失患者，显然我们将严重低估这群病患的平均生存时间。

from lifelines.plotting import plot_lifetimes
from numpy.random import uniform, exponential

N = 25

CURRENT_TIME = 10

actual_lifetimes = np.array([
    exponential(12) if (uniform() < 0.5) else exponential(2) for i in range(N)
])
observed_lifetimes = np.minimum(actual_lifetimes, CURRENT_TIME)
death_observed = actual_lifetimes < CURRENT_TIME

ax = plot_lifetimes(observed_lifetimes, event_observed=death_observed)

ax.set_xlim(0, 25)
ax.vlines(10, 0, 30, lw=2, linestyles='--')
ax.set_xlabel("time")
ax.set_title("Births and deaths of our population, at $t=10$")
print("Observed lifetimes at time %d:\n" % (CURRENT_TIME), observed_lifetimes)

Observed lifetimes at time 10:
[10.   1.1   8.   10.  3.43   0.63   6.28   1.03   2.37   6.17  10.
   0.21   2.71   1.25  10.   3.4  0.62   1.94   0.22   7.43   6.16  10.
   9.41  10.  10.]

此外，如果简单地对在第10个月观察到的患者（包括右删失患者）生存时间取平均值，也会低估正式的生存时间。下图我们绘制了这群病患的真实生存时间。

ax = plot_lifetimes(actual_lifetimes, event_observed=death_observed)
ax.vlines(10, 0, 30, lw=2, linestyles='--')
ax.set_xlim(0, 25)

通过以上距离，我们可以认识到了右删失数据对生存时间估计的影响很大。生存分析方法就是用来解决这个问题的。当然啦，当不存在任何删失样本的情况，生存分析方法也是适用且很好用的。

上述例子，为了便于理解，我们假设在t=0时刻开始观察。实际上，观察时间并不总是在0时刻开始。例如这种情境：在一个小卖部中，t=0时刻定义为小卖部早上开门营业那一时刻，“出生”事件定义为顾客走进小卖部——顾客可以在任意时刻走进小卖部，不一定是0时刻。在生存分析中，生存时间是相对的，因为样本可能在不同时刻开始。（实际上，在生存分析中，我们对样本开始时刻、结束时刻都不感兴趣，而是对从开始到结束的时间长度感兴趣）

接下来介绍生存分析中的两个基础概念：生存函数和风险函数。

生存函数

记T为时间，例如，一对夫妻结婚多久了、再如用户进入一个网页需花费的时间（如果网页404了，则该时间为无限长）。生存函数S（t）定义为：
$$S(t)=Pr(T>t)$$
通俗地说，生存函数定义了“死亡”事件直到t时刻都还没有发生的概率，或者说，生存时间大于t的概率。
生存函数具有以下性质：

1. $0\le S(t)\le 1$
2. $F_T(t)=1-S(t)$,其中$F_T(t)$是$T$的CDF
3. $S(t)$是$t$的非增函数。

风险曲线

除了生存时间，我们还对已知“死亡”事件在$T$时刻前未发生的前提下，在$t$时刻发生的概率感兴趣。将该概率表示为数学表达式，即：
$$\underset{\delta \to 0}{\lim }Pr(t\le T\le t+\delta t|T>t)$$
当$\delta t$趋于0的时候，$Pr(t\le T\le t+\delta t|T>t)$将趋于0，我们将无法考察概率，因此我们将上式除以$\delta t$,这就是$t$时刻的风险函数$h(t)$的定义:
$$h(t)=\underset{\delta t\to 0}{\lim }\frac{Pr(t\le T\le t+\delta t|T>t)}{\delta t}$$
不难推导，得到：
$$h(t)=\frac{-S^{-1}(t)}{S(t)}$$
酷！这是一个很容易求解的常微分方程！解该方程，得到：
$$S(t)=exp(-{\int }_0^{t}h(z)dz)$$
上式的有趣之处在于，它定义了生存函数。有了这个关系式，可以发现$S(t)$和$h(t)$之间具有唯一的关系，两者可以相互转换。
然后，可以得到，“死亡”事件在$t$时刻发生的概率密度如下：
$$f_T(t)=h(t)exp(-{\int }_0^{t}h(z)dz)=h(t)S(t)$$

当然，我们无法观察到真正的生存曲线，只能通过观察到的数据去估计。估计生存函数和风险比的方法有很多，详情参考这个链接：用lifelines来做预测（尚未上线）。