C#

使用生存分析

Zvi Topol

生存分析 (SA) 是一门学科的统计数据,重点是估计到的事件的时间。 你通常会适用于临床研究,以帮助确定某些药物 (时间到病人死亡) 的有效性、 可靠性的软件系统 (失效时间) 和信贷分析 (时间到拖欠贷款) 的生存分析方法。

涉及两组病人的医药临床研究是如何工作的最好的例子。 对照组成员管理的安慰剂。 测试组成员收到针对这种疾病的实验性药物。 然后应用生存分析方法来确定是否有统计学意义在病人的存活,两组间差异。 "时间到"事件在这种情况下是从研究开始直到病人死亡的时间。

若要公开你使用 SA,我掩护和 C# 实现的一种常用的估计方法,称为卡普兰 Meier 估计的基本概念。 我将使用估算的移动应用程序的生存概率的现实世界的例子。

想象一个软件开发公司产生两个单独移动应用程序,名为 X 和 Y。 每个单独的团队被开发。 该公司是渴望学习如何强大的移动应用程序,并确定是否一个应用程序是显著较不可靠,需要更多努力,提高其可靠性。

在任何给定的时刻,可以有很多实例的 X 和 Y 还活着并在客户的移动设备上运行。 因此,移动应用程序崩溃是最有趣的是。 在这种情况下,长时间对这一事件,将表明应用程序更为健壮或有生存力更强。

在演示程序中,我第一次将显示 X 和 Y 的移动应用程序的生存数据 (请参见图 1)。 数据显示两个应用程序的 10 个不同用户用 Id 从零到九。 在我的示例中,应用程序可以要么崩溃 (由事件描述 = 在截图中的应用程序崩溃) 或获取由用户关闭 (由事件描述 = 关闭的应用程序)。 此外可以记录的事件发生的那一天。

生存分析演示显示移动应用程序的生命周期
图 1 显示的移动应用程序生命周期的生存分析演示

SA 的基本概念

在 SA 中最基本的概念是生存函数。 这通常是由 s (t) 表示。 如果 X 是一个随机变量 (变量的值是基于机会的成果) 表示的事件的时间,然后 s (t) = Pr (X > t)。 换句话说,s (t) 是时间 t 后生存的可能性。 套装 (0) 被定义为 1。 生存函数被有关使用寿命分布函数。 这通常由 f (t) 表示,定义为 f (t) = Pr (X<= t),换句话说,直到时间 t 发生的事件的概率。 因此,f (t) = 1 — — s (t)。 然后定义事件密度函数 f (t) 是 dF (t) / dt — — f (t),如果 f (t) 是可微的第一阶导数。 因此,f (t) 可以作为率 (单位时间) 事件的思路。

风险函数,另一个基本概念,等于 f(t)/S(t),是针对那些还活着,时间 t 的个人,时间 t 的事件率。

您可以指定生存函数参数,使用显式函数族。 你也可以推断他们非参数地从现有的数据,而无需关闭的参数化窗体。 半参数的规范,是参数和非参数规范之间的组合,也是可能的。 指数分布是简单和受欢迎的参数化功能家庭描述生存功能由于其具有吸引力的数学属性。

例如,s (t) = exp(-0.05t) 是一个生存在中绘制实时指数分布从图 2。 生存函数的形式 s (t) = exp(-at) (哪里是控制造成的危害率参数可以描述这种分配)。 分布函数 f (t) 由给的生存期 = 1 — — s (t) = 1 — — exp(-at)。 图 2 帮助我们直观显示随着时间的推移生存函数的行为方式。

随着时间的推移生存功能的工作方式
随着时间的推移生存功能的工作方式的图 2

使用给定的参数化模型,您可以使用实际数据来估计模型的参数。 指数分布,则该参数。 这样做的一种方法是使用最大似然估计 (MLE) 方法,但这完全是另一个主题。

我将重点实施非-­生存函数的参数估计。 这就是,我不会设置为生存函数的预定义的模型和估计模型参数。 相反,我会直接从观察到的数据派生的生存函数。 描述了如何这样做之前,我要解释的 SA 叫截尾的另一个重要概念。

截尾时发生一些意见集中的数据并不完整。 有些时候,你都忘了项目观察。 在我的示例中,这将意味着移动应用程序结束它的执行不会导致崩溃 (投掷致命的异常)。 应用程序已由用户正常关闭。 尽管可以有其他应用程序结束时,没有崩溃的原因,我会认为应用程序崩溃或获取由用户关闭。

有两个主要风味的审查 — — 右截尾和左截尾。 右截尾时发生的开始时间已知的但事件的时间丢失。 左截断数据时发生的事件时间存在,但开始时间是失踪。 右截尾是出现在我的示例。

使用的卡普兰 Meier 估计估计生存函数

卡普兰 Meier (公里) 估计是一种非参数化的算法,估计生存函数。 派生的 KM 估计涉及的高等数学,包括鞅理论和计票过程,使用和超出了本文的范围。 实施 KM 估计,但是,是简单的和基于计数。

考虑计算 KM 估计的移动应用程序 X 的生存。 KM 估计需要跟踪的三个不同的计数:

  1. 移动应用程序 X 的多少个实例都仍然运行。 这是在我执行使用变量 atRisk 来表示。
  2. 坠毁的实例数。 这被跟踪在坠毁的变量中。
  3. 优雅地完成执行的实例数。 这些计数使用送检的变量。

下面的代码行 (移动应用程序 X) 正在使用 CrashMetaData 类进行编码中表示的生存数据图 3

var appX = new CrashMetaData[] {new CrashMetaData{UserID = 0,
  CrashTime = 1, Crashed = false},
           new CrashMetaData{UserID = 1, 
              CrashTime = 5, Crashed = true}, 
           new CrashMetaData{UserID = 2, 
              CrashTime = 5, Crashed = false},
           new CrashMetaData{UserID = 3, 
              CrashTime = 8, Crashed = false},
           new CrashMetaData{UserID = 4, 
              CrashTime = 10, Crashed = false},
           new CrashMetaData{UserID = 5, 
              CrashTime = 12, Crashed = true},
           new CrashMetaData{UserID = 6, 
              CrashTime = 15, Crashed = false},
           new CrashMetaData{UserID = 7, 
              CrashTime = 18, Crashed = true},
           new CrashMetaData{UserID = 8, 
              CrashTime = 21, Crashed = false},
           new CrashMetaData{UserID = 9, 
              CrashTime = 22, Crashed = true}};

图 3 生存数据的移动应用程序 X

用户 Id 坠毁 检察
0 1   X
1 5 X  
2 5   X
3 8   X
4 10   X
5 12 X  
6 15   X
7 18 X  
8 21   X
9 22 X  

生存的数据包含在几天 (由 CrashTime 编码) 和有关事件是否指应用程序的信息的事件时间崩溃或审查。 如果 Crashed 等于 true,应用程序崩溃。 否则,应用程序正确关闭 (换句话说,进出都要审查)。 此外,用户 Id 的字段跟踪应用程序的实例。

KM 估计是在 EstimateKaplanMeier 方法中实现的。 此分区数据的不同非重叠的时间间隔基于时间段 (在我的示例中这是应用程序崩溃) 的事件。 它跟踪的每个时间间隔中的计数。

它是重要的是要注意数到多少应用程序仍上升了,只是之前的事件 (这是由于计票过程的数学公式) 进行运行。 所以在我的示例,其中包括 0 到 5 天,第一次间隔中 9 出的 10 个实例是向上和运行只是前一天 5 (完成 1 次运行一个实例)。 在直至并包括第 5 天的时间间隔,我有一个崩溃 (其中定义的时间间隔) 和 2 实例整理 (在 1 和 5 天)。 参见图 4

一天时间间隔由 KM 估计创建
图 4 天的时间间隔创建的 KM 估计

生存函数的 KM 估计数然后是生存的产品在分区中的计数从派生的所有不同的时间间隔:

1 — — (间隔中崩溃)/(那些在风险只是在该间隔结束之前)

EstimateKapalanMeier 方法返回一个 SurvivalCurve 类的对象。 这表示估计的生存函数。 输出是单步执行函数。 每个步骤是生存函数的值在相应时间间隔 (如估计 KM 估计)。 图 5 包括生存分析演示程序输出对应的 SurvivalCurve 对象 (的 X 和 Y 两个应用程序) 的一部分。

生存分析演示输出为 KM 估计为 X 和 Y 的应用程序的
图 5 生存分析演示输出为 KM 估计为 X 和 Y 的应用程序的

图 6 包括估计为 X 的移动应用程序的步骤生存函数的一个阴谋。 在剧情中的每一步的短竖线表示在对应于该步骤的间隔期间崩溃事件的多个匹配项。

KM 的函数估计的生存为 X 的移动应用程序
图 6 公里的函数估计的生存为 X 的移动应用程序

然后可以使用估计值来推断的中位生存时间或由该实例的一半将活着的时间。 这应该发生在一些点在 12 天之间的时间 (生存概率估计在哪里 0.711 > 0.5) 和 18 (生存概率在哪里 0.474 < 0.5)。 有几种方法在 SA 文学描述如何确切地计算此数量,因为它通常落在两个步骤之间。

我会将中位生存时间定义为最小存活时间为其生存函数小于 0.5,导致为 X 的移动应用程序中中, 位生存时间为 18 天。 这个数量的解释是,由天 18 的移动实例的应用程序 X 崩溃和一半留一半就和运行。 此实现计算中位生存时间在 GetMedianSurvivalTime 方法中。

你可以回答使用 KM 估计数的另一个问题是是否有生存能力的两个 (或多个) 不同的应用程序的差异。 处理这一问题的一种方法是以可视方式绘制的 KM 估计对应于每个应用程序。 这种类型的情节描述在图 7,并将应用程序的 X 和 Y 的估计的生存功能进行比较。

KM 估计为 X 和 Y 的移动应用程序的
图 7 公里估计为 X 和 Y 的移动应用程序

绿色曲线表示生存函数的应用程序 X,蓝色曲线表示生存函数的应用程序 Y。

从情节,你可以看到生存功能的应用程序 X 上衣生存功能的应用程序 Y。 因此,你可以推断应用程序 X 有比应用程序 Y 生存力更强,因此,是更强健。

虽然可视化生存函数可能有助于确定可生存性的差异,某些情况下不是一样很明确的。 幸运的是,有统计的方法,这种差异在正规、 严格的方式,称为日志等级测试中的测试。 这是一种算法,测试是否有显著差异的非参数化方式的两个 (或更多) 生存分布。 SA 文学包括这个问题的详细的讨论,而大多数 SA 统计库包括日志等级测试实现。

值得一提的还有另一种流行算法在非参数化方式叫纳尔逊 · 阿隆 (NA) 估计估计生存函数。 NA 估计从生存数据的累积风险函数。 然后,您可以从这个估计使用一个数学公式,它关系到累积风险函数派生生存函数。 SA 文献中,可以找到有关此估计的更多详细信息。

总结

我从生存分析的统计分支介绍基本概念和术语。 我展示了如何执行非参数化的卡普兰 Meier 估计,并将其应用到一个比较稳健的移动应用程序的示例。 这个估计可以帮助确定是否有生存能力的两个应用程序的差异。 我也提到严格的统计测试,以检查存在差异,称为日志等级测试。 另一个数量的 KM 估计推算是中位生存时间,也指向应用程序 X 和 Y 之间的可生存性差异。

最后,我提到的纳尔逊 Aalen 估计为替代性的非参数方法估算的生存函数。 它并不直接估计生存函数像 KM 估计,但是,而是估计累积风险函数。 然后,您可以从累积风险函数派生的生存函数。

这只触及表面的丰富领域的 SA。 应用程序跨领域,从医学到在许多统计软件包中实现工程和其方法和算法。 与移动应用程序和软件作为服务企业部署扩散,我预计 SA 方法可以在发挥作用监测和提高质量的这种部署。

Zvi Topol  工程作为一个在市场营销分析在纽约城的高级科学家。他设计和适用非线形大规模优化算法和统计方法,以改善大财富 500 强企业的市场营销规划。

衷心感谢以下技术专家对本文的审阅:博士。 詹姆斯 · 麦卡弗里 (Microsoft 研究)