数据科学工程实践:用户行为分析与建模、A/B实验、SQLFlow
上QQ阅读APP看书,第一时间看更新

2.2 生存分析的理论框架

生存分析包括4个主要过程,如图2-2所示,本节将逐一进行介绍。

000

图2-2 生存分析的4个主要过程

2.2.1 生存分析基本概念界定

1. 事件

事件分为起始事件和终点事件。起始事件即生存分析的起点事件,所有研究对象的起始事件均相同;终点事件指研究关注的具体事件,生存分析中部分研究对象可以观察到关注事件发生,能够获取准确的发生时间,为分析提供完全信息,本章后续讨论中提到的“事件”均指“终点事件”。

在二手车定价案例中,起始事件是平台发布车辆出售信息,终点事件是用户下单购买二手车。

2. 生存时间

生存时间不关注事件发生的客观时间,只关注从起始事件开始到终点事件发生之间的时间间隔。例如,病人治愈所花费的时间;首婚人群的婚姻持续时间;某个系统在故障前良好运转的时间等。生存时间不呈正态分布,因此不能用生态分布假设对生存时间分布参数进行估计。

在二手车定价案例中,车辆生存时间是从平台发布信息到车辆被售出的时间间隔,也代表了车辆在库存存放的时长。

3. 删失

某些研究对象在观察时间窗内无法获取事件发生时间,我们将这种情况称为删失,其中以右删失最为常见,右删失是指已知研究对象的观察起始时间,但无法获取终点事件发生的具体时间,这一现象是由以下几个原因导致的。

  • 研究对象在观察时间窗内还未发生有效事件。
  • 研究对象在观察时间窗内由于某些原因被丢失。
  • 研究对象在事件发生前由于非事件原因脱离有效观测。

在二手车定价案例中,因调价导致原始价格下车辆后续销售状况不明就是一种右删失,其原因可以归纳为第三类。二手车案例数据删失举例如图2-3所示。

000

图2-3 二手车案例中数据删失举例

针对车辆3,有效的处理方式是将数据拆分为调价前与调价后两段,调价前时间段由于受调价影响,无法追踪车辆最终是否出售,因此定义为删失,调价后在新价格阶段能够有效追踪销售情况,具体说明如图2-4所示。

000

图2-4 二手车案例中调价信息处理

除右删失外,还有左删失与区间删失。左删失是指确定研究对象在某一时刻之前发生了有效事件,但发生具体时间不详;区间删失是指已知某一研究对象在某一时间段内发生了有效事件,但发生具体时间不详。

4. 风险中数量

风险中数量指在观察时间窗内可追踪其状态且未发生事件的对象数量,但不包括如下两项。

  • 截止到当前时间已经发生了事件/流失的对象。
  • 截止到当前时间已经右删失的对象(由于非事件原因脱离观测的对象)。

在二手车定价案例中,假设在初始价格下,100辆车在信息发布后的20天内,有40辆车被售出,有4辆车在售出前进行了调价,那么在第20天末,风险中数量是100-40-4=56个。

5. 生存函数

对于观察时间窗内的任意时刻tt>0),生存函数反映的是研究对象到该时刻仍未发生事件的概率。生存函数是每个时刻生存概率的乘积,故也称为累积生存概率函数。

生存函数的公式表示为

S(t) = P(T > t) = 1 − F(t)

F(t)代表生存时间的累积分布函数,表示事件发生时间未超过时刻t的概率。

生存函数具备以下3个特性。

  • S(t)∈[0,1],且S(t)单调递减。
  • 在起始时刻t=0时,所有对象均处于存活状态,此时S(t)=1。
  • t趋于无穷大时(t=∞),生存概率趋近于0,S(t)=S(∞)=0。

当生存时间t为连续型随机变量时,生存函数表示为

000

其中,f(t)为概率密度函数,是F(t)的导数。

基于生存函数可以绘制生存曲线。生存函数对应一条从1到0下降的曲线。曲线越靠左越陡峭,代表生存率越低或生存时间越短,一般的生存曲线如图2-5所示。

000

图2-5 生存曲线示意图(非案例数据)

在二手车定价案例中,利润函数中的P(t,p)即为生存函数,代表在价格等于p的条件下,车辆存放到第t日仍然没有卖出去的概率。每个时刻生存概率越大,代表车辆出售的速度越慢。基于P(t,p)可绘制出不同价格条件下的生存曲线,而000对应的是曲线下的面积,即车辆的平均出售时长。我们的目标就是通过绘制不同价格、不同参数条件下车辆的生存曲线,求解车辆平均出售时长,进而制定出最优价格策略。

6. 风险函数与累积风险函数

风险函数也可以称为条件死亡率,指的是在时间t之前未发生任何事件而恰好在时间t发生事件的概率。

风险函数的公式表示为

000

当生存时间t为连续型随机变量时,风险函数表示为

000

累积风险函数公式为

000

生存函数与风险函数的关系表示为

000

7. 半衰期

半衰期也可以称为中位生存时间,指恰好一半个体未发生终点事件的时间。由于删失的存在,无法直接以平均时间反应事件发生的时间水平,因此生存分析在描述生存时间时一般以半衰期,也就是存活概率50%时对应消耗的时间来表述。同理,也可以基于需要计算存活概率到其他百分位数量的消耗时间,半衰期示例如图2-6所示。

000

图2-6 半衰期(中位生存时间,非案例数据)

在二手车定价案例中,半衰期指恰好还有50%的车辆没有卖出去对应的时间,由于生存概率是一个随时间变化的值,无法直接进行对比,因此通常可以采用半衰期作为反映生存概率高低的基础指标,通过对半衰期的监控及对比反映销售情况。

2.2.2 生存函数刻画及简单对比

2.2.1节介绍了生存函数与风险函数的定义及公式,本节将介绍对两个函数进行刻画的具体方法。目前,对生存函数及风险函数的刻画方法分为非参数和参数两类。其中,最常用的是非参数方法中对生存函数进行刻画的KM曲线法和对累积风险函数进行刻画的Nelson-Aalen曲线法。本小节会着重介绍KM曲线法,然后对其他方法做简要说明。

1. KM曲线法

KM曲线法作为一种非参数方法,不对数据分布做任何假设,而是直接用概率乘法定理估计生存率。这一方法的优势在于能够直观地观察生存曲线,便于不同生存曲线之间进行简单对比,但无法建立数学模型对多个影响因素进行分析。

假设我们有(t1<t2<…<tk)共k个观测时刻及N个样本,dj代表在tj时刻发生事件的人数,mj表示在(tj,tj+1)时间段内删失的人数,那么tj时刻依然处于风险中的人数可以表示为

nj = (mj + dj) + … + (mk + dk)

tj时刻的风险率为dj/nj,KM曲线对应生存函数的估计表示为

000

KM曲线的估计可以用经典Greenwood公式进行计算,从而得到生存函数的置信区间,具体计算公式为

000

2. KM曲线对比

我们可以通过对研究对象进行分组,绘制多条生存曲线,但只通过直接观察无法确定曲线之间是否具有显著性差异,还需要引入严谨的统计检验方法。其中,对数秩检验在生存曲线的比较中应用最为广泛。

对数秩检验是一种非参数检验方法,检验的原假设是不同分组之间的生存率没有显著差异。在原假设为真的条件下,对数秩检验以整体的风险概率作为理论风险概率计算理论事件数,通过将真实事件数与理论事件数进行对比,判断生存曲线之间是否存在显著差异。以两个分组下的生存曲线对比为例,对数秩检验的具体步骤如下。

  • 将两组数据混合后按照生存时间排序,通过KM曲线法计算合并后数据的整体风险率及生存率。
  • 以整体风险率作为理论风险率,计算每组数据在原假设成立的情况下,各时刻的理论事件数(期望事件数)。
  • 对各组的理论事件数进行求和,并计算检验统计量量,在原假设成立时,统计量服从自由度为组数减1的卡方分布,具体检验公式表示为
  • 000
  • 依据P值判断原假设是否成立,p-value≤0.05代表不同组生存曲线之间存在显著差异,反之则没有显著差异。

由此可见,对数秩检验实际上是一种单因素分析,该方法能够有效实现单因素分类后的组间比较,但无法有效实现多因素分析。

3. Nelson-Aalen累积风险曲线

与KM曲线相同,Nelson-Aalen累积风险曲线法也是一种非参数方法。该方法适用于对删失数据的处理,曲线基于观测数据对风险函数及累积风险函数进行刻画。

其中,累积风险函数的估计表示为

000

4. 生存函数的参数估计

参数方法假定生存时间符合某种分布,根据样本观测值来估计假定分布模型中的参数,以获得生存时间的概率密度模型。通常,假定生存时间服从的分布主要有指数分布、威布尔分布、对数正态分布及对数logistic分布等,下面对指数分布及威布尔分布做简要介绍。

指数分布的生存函数表示为

S(t) = e-λt

风险函数表示为

h(t) = λ

可见,指数分布假设时事件发生时间完全随机,风险函数与时间无关。

威布尔分布的生存函数表示为

S(t) = e(-λtα)λ>0,α>0

风险函数表示为

h(t) = λα(λt)α-1

其中,λ为尺度参数,决定分布的分散度;α为形状参数,决定分布的形态。α=1时为指数分布,威布尔分布允许风险函数单调递增或递减。

2.2.3 生存函数回归及个体生存概率的预测

1. Cox比例风险回归模型

前文提到,KM曲线之间的对比只能在单一分类变量之间进行,并且KM曲线只描述了该单变量和生存时间之间的关系而忽略了其他变量的影响。为了解决这个问题,Cox比例风险回归模型诞生了。Cox比例风险回归模型通过对风险率进行估计,实现了对包含分类变量及连续变量的多变量回归,并控制了其他变量下,单个变量对生存概率的影响。此外,基于回归的结果可以实现在给定条件下个体生存概率的预测。

Cox比例风险模型公式为

h(t;x)=h0(t)eβx

与其他回归分析相同,生存回归模型可以包括多个自变量,假设有一个影响风险率的自变量,则公式可具体表述为

h(t;x1,x2,…xp)=h0(t)eβ1x1+β2x2+…+βpxp

对于风险函数而言,基准风险率h0(t)是所有自变量取值均为0的风险率,(β1,β2,…,βp)是模型的偏回归系数,是一组带估计的参数,反应了自变量对风险率的影响。模型基于恒定比例风险假设,在不同时刻,相同水平下风险率相对基准风险率的比例固定,不随时间变化,模型采用指数函数形式可以确保风险率大于0。与一般的回归分析不同,Cox模型不是直接用生存时间作为回归方程的因变量,自变量对生存时间的影响通过风险函数与基准函数的比值来表示,其中模型中h0(t)为非参数部分,β为需要估计的参数,因此Cox比例风险回归模型是一种半参数模型。

我们对模型进行解读:当自变量xp为连续型变量时,偏回归系数βp可表示为,在除xp之外其他变量均不变的情况下,xp每增加一个单位,风险率变化eβp倍;当自变量xp为离散型变量时,假设xp有0和1两个分类,则偏回归系数βp可表示为,在除xp之外其他变量均不变的情况下,类别1相对于类别0的风险率是eβp倍。

总体来说,对于影响因子的解读可以分为以下3类。

  • β > 0时,风险率随着x的增加而增加,我们称之为危险因素。
  • β < 0时,风险率随着x的增加而减少,我们称之为保护因素。
  • β = 0时,风险率不随x的变化而变化,我们称之为无关因素。

2. Cox比例风险回归模型的估计

Cox比例风险回归模型中的偏回归系数采用部分似然函数进行估计。假设我们有一个研究对象,生存时间按照由大到小排序:

t1t2 ≤ … ≤ tn

将在时刻ti依然处于风险中的所有对象构成危险集,记为R(ti)。随着ti的增大,危险集里的数量逐渐减少直到消失,部分似然函数表示为

000

对于研究对象i,其在ti时刻发生事件的概率是两部分的乘积:一部分是在ti时刻之前依然处于风险中的概率;另一部分是在ti时刻的风险集中,恰好第i个研究对象发生事件的概率。PL(β)没有考虑前者,因此成为部分似然函数。假设对象iti时刻发生事件,则记δi=1,否则δi=0,对PL(β)取对数得到:

000

求关于βj(j=1,2,…,p)的一阶偏导,并求其等于0的解,即可得到βj的最大似然估计值。

在估计出模型中的β值之后,将每个个体的自变量取值带入公式,就可以得到每个个体的生存函数取值,绘制个体的生存函数曲线,从而实现对个体生存概率的预测。