生存分析的名称源于临床研究，其中预测死亡时间，即生存，通常是主要目标。

生存分析是一种回归问题（人们想要预测一个连续值），但有一个转折点。它与传统回归的不同之处在于，在生存分析中，结果变量既有一个事件，也有一个与之相关的时间值，部分训练数据只能被部分观察——它们是被删失的。本文用R语言生存分析晚期肺癌患者数据（查看文末了解数据获取方式）。

普通最小二乘回归方法不足，因为事件发生的时间通常不是正态分布的，并且模型无法处理删失，但这在生存数据中很常见。

为什么要做生存分析：右删失

在某些情况下，可能无法观察到事件时间：这通常称为 右删失。在以死亡为事件的临床试验中，当发生以下情况之一时，就会发生这种情况。1。当一定数量的参与者死亡时，研究结束。2。参与者退出研究。3。 研究达到预定的结束时间，并且一些参与者存活到结束。在每种情况下，幸存的参与者离开研究后，我们都不知道他们会发生什么。然后我们有一个问题：

当对于某些个体，我们只观察到他们的事件时间的下限时，我们如何对经验分布进行建模或进行非负回归？

上图说明了右删失。对于参与者 1，我们看到他们何时死亡。参与者 2 退出了，我们知道他们一直活到那时，但不知道后来发生了什么。对于参与者 3，我们知道他们活到了预定的研究结束，但又不知道之后发生了什么。

生存函数和风险函数

生存分析中的两个关键工具是生存函数和风险函数。

生存函数：它是一个函数，用于给出我们有兴趣知道的任何对象是否会在任何指定时间之后存活的概率。在数学上它可以由以下公式表示

其中 S(t) 是一个生存函数，其中 T 是一个连续随机变量，是一个事件的时间。F(t) 是区间[0,∞) 上的累积分布函数。

我们也可以用风险函数来写生存函数。假设事件尚未发生 ，风险率λ(t) 是事件在时间t发生的瞬时概率的主要值。

那么关键问题是如何估计风险和/或生存函数。

Kaplan Meier的非参数估计

在非参数生存分析中，我们要估计生存函数没有协变量，并且有删失。如果我们没有删失，我们可以从经验 CDF 开始. 这个等式简洁地表示：

有多少人随着时间的推移而死亡? 那么生存函数就是：还有多少人还活着？

但是，我们无法回答一些人被时间t删失时提出的这个问题.

虽然我们不一定知道有多少人在任意时间t幸存下来，我们知道研究中有多少人仍然处于风险之中。我们可以使用它来代替。将学习时间划分区间, 其中每个ti是参与者的事件时间或删失时间。假设参与者只能在观察到的事件时间失效。假设没有人在同一时间死去（没有关系），我们可以查看每次有人死去的时间。我们说在那个特定时间死亡的概率是，并说在任何其他时间死亡的概率是0.

在温和的假设下，包括参与者具有独立且相同分布的事件时间，并且删失和事件时间是独立的，这给出了一个一致的估计量。上图给出了一个简单案例的 Kaplan Meier 估计示例。

生存分析用于各种领域

例如：

• 用于患者生存时间分析的癌症研究，
  
• “事件历史分析”的社会学，
  
• 在工程中用于“故障时间分析”。
  

在癌症研究中，典型的研究问题如下：

• 某些临床特征对患者生存有何影响
  
• 一个人能活3年的概率是多少？
  
• 患者组之间的生存率是否存在差异？
  

第1部分：生存分析简介

本演示文稿将介绍生存分析 ，参考：

Clark, T., Bradburn, M., Love, S., & Altman, D. (2003). Survival analysis part I: Basic concepts and first analyses. 232-238. ISSN 0007-0920.

我们今天将使用的一些软件包包括：

• lubridate
  

library(survival)

什么是生存数据？

事件时间数据由不同的开始时间和结束时间组成。

癌症的例子

• 从手术到死亡的时间
  
• 从治疗开始到进展的时间
  
• 从响应到复发的时间
  

其他领域的例子

事件发生时间数据在许多领域都很常见，包括但不限于

• 从艾滋病毒感染到艾滋病发展的时间
  
• 心脏病发作的时间
  
• 药物滥用发生的时间
  
• 机器故障时间
  

生存分析别名

由于生存分析在许多其他领域很常见，因此也有其他名称

• 可靠性分析
  
• 持续时间分析
  
• 事件历史分析
  
• 事件发生时间分析
  

肺数据集

数据包含来自北中部癌症治疗组的晚期肺癌患者。今天我们将用来演示方法的一些变量包括

• 时间：以天为单位的生存时间
  
• 状态：删失状态1 =删失，2 =失效
  
• 性别：男= 1女= 2
  

删失类型

某个主题可能由于以下原因而被删失：

• 后续损失
  
• 退出研究
  
• 固定学习期结束前没有活动
  

具体来说，这些是删失的示例。

分配随访时间

• 受删失的主题仍会提供信息，因此必须适当地包含在分析中
  
• 随访时间的分布存在偏差，在接受检查的患者和有事件的患者之间可能有所不同
  

生存数据的组成部分

对于主题ii：

• 活动时间Ti
  
• 删失时间Ci
  
• 事件指标δi：
  
• 1，如果观察到的事件（即  Ti≤CiTi≤Ci）
  
• 如果检查，则为0（即  Ti>CiTi>Ci）
  

9. 观测时间Yi=min(Ti,Ci)Yi=min(Ti,Ci)
   

lung数据中提供了观察时间和事件指示

• 时间：以天为单位的生存时间（YiYi）
  
• 状态：删失状态1 =删失，2 =死亡（δiδi）
  

在R中处理日期

数据通常带有开始日期和结束日期，而不是预先计算的生存时间。第一步是确保将这些格式设置为R中的日期。

让我们创建一个小的示例数据集，其中sx_date包含手术日期和last_fup_date上次随访日期的变量。

date_ex <- 
  tibble(
    sx_date = c("2007-06-22", "2004-02-13", "2010-10-27"), 
    last\_fup\_date = c("2017-04-15", "2018-07-04", "2016-10-31")
    )
date_ex

## # A tibble: 3 x 2
##   sx\_date    last\_fup_date
##   <chr>      <chr>        
## 1 2007-06-22 2017-04-15   
## 2 2004-02-13 2018-07-04   
## 3 2010-10-27 2016-10-31

我们看到它们都是字符变量，通常都是这种情况，但是我们需要将它们格式化为日期。

格式化日期-基数R

date_ex %>% 
  mutate(
    sx\_date = as.Date(sx\_date, format = "%Y-%m-%d"), 
    last\_fup\_date = as.Date(last\_fup\_date, format = "%Y-%m-%d") 
    )

## # A tibble: 3 x 2
##   sx\_date    last\_fup_date
##   <date>     <date>       
## 1 2007-06-22 2017-04-15   
## 2 2004-02-13 2018-07-04   
## 3 2010-10-27 2016-10-31

• 请注意，R格式必须包含分隔符和符号。例如，如果您的日期格式为m / d / Y，则需要format = "%m/%d/%Y"
  

格式化日期-lubridate程序包

我们还可以使用该lubridate包来格式化日期。在这种情况下，请使用ymd功能

date_ex %>% 
  mutate(
    sx\_date = ymd(sx\_date), 
    last\_fup\_date = ymd(last\_fup\_date)
    )

## # A tibble: 3 x 2
##   sx\_date    last\_fup_date
##   <date>     <date>       
## 1 2007-06-22 2017-04-15   
## 2 2004-02-13 2018-07-04   
## 3 2010-10-27 2016-10-31

• 请注意，与基本R选项不同，不需要指定分隔符
  

计算生存时间

现在日期已格式化，我们需要以某些单位（通常是几个月或几年）计算开始时间和结束时间之间的差。在base中R，用于difftime计算两个日期之间的天数，然后使用将其转换为数字值as.numeric。然后将除以365.25年的平均天数转换为年。

date_ex %>% 
  mutate(
    os_yrs = 
      as.numeric(
        difftime(last\_fup\_date, 
                 sx_date, 
                 units = "days")) / 365.25
    )

## # A tibble: 3 x 3
##   sx\_date    last\_fup\_date os\_yrs
##   <date>     <date>         <dbl>
## 1 2007-06-22 2017-04-15      9.82
## 2 2004-02-13 2018-07-04     14.4 
## 3 2010-10-27 2016-10-31      6.01

计算生存时间

操作员可以%--%指定一个时间间隔，然后使用将该时间间隔转换为经过的秒数as.duration，最后除以dyears(1)，将其转换为年数，从而得出一年中的秒数。

## # A tibble: 3 x 3
##   sx\_date    last\_fup\_date os\_yrs
##   <date>     <date>         <dbl>
## 1 2007-06-22 2017-04-15      9.82
## 2 2004-02-13 2018-07-04     14.4 
## 3 2010-10-27 2016-10-31      6.02

事件标标

对于生存数据的组成部分，我提到了事件指示器：

事件指标δiδi：

• 1，如果观察到的事件（即  Ti≤CiTi≤Ci）
  
• 如果检查，则为0（即  Ti>CiTi>Ci）
  

在lung数据中，我们有：

• 状态：删失状态1 =删失，2 =失效
  

生存函数

受试者可以存活超过指定时间的概率

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

S(t)S(t)：生存函数F(t)=Pr(T≤t)F(t)=Pr(T≤t)：累积分布函数

理论上，生存函数是平滑的；在实践中，我们以离散的时间尺度观察事件。

生存概率

• 生存概率在某个时间，S(t)S(t)，是存活超过该时间，考虑到个体已存活刚刚在此之前，时间的条件概率。
  
• 可以估计为当时活着但没有损失的随访患者人数除以当时的活着患者人数
  
• 生存概率的Kaplan-Meier估计是这些条件概率的乘积
  
• 在时间0，生存概率为1，即  S(t0)=1S(t0)=1
  

创建生存对象

Kaplan-Meier方法是估计生存时间和概率的最常用方法。这是一种非参数方法，可产生阶跃函数，每次事件发生时，阶跃下降。

• 创建一个生存对象。对于每个主题，将有一个条目作为生存时间，+如果主题是经过删失的，则后面跟一个。让我们看一下前10个观察值：
  

##  \[1\]  306   455  1010+  210   883  1022+  310   361   218   166

用Kaplan-Meier方法估算生存曲线

• 该survfit函数根据公式创建生存曲线。让我们为整个同类群组生成总体生存曲线，将其分配给object f1，然后查看names该对象的：
  

names(f1)

##  \[1\] "n"          "time"       "n.risk"     "n.event"    "n.censor"  
##  \[6\] "surv"       "std.err"    "cumhaz"     "std.chaz"   "start.time"
## \[11\] "type"       "logse"      "conf.int"   "conf.type"  "lower"     
## \[16\] "upper"      "call"

该survfit对象将用于创建生存曲线的一些关键组件包括：

• time，其中包含每个时间间隔的起点和终点
  
• surv，其中包含每个对应的生存概率 time
  

Kaplan-Meier图

现在， 绘制对象 获得Kaplan-Meier图。

plot(survfit(Surv(time, status) ~ 1, data = lung),

• 基数R中的默认图显示了具有相关置信区间（虚线）的阶跃函数（实线）
  
• 水平线代表间隔的生存时间
  
• 时间间隔由事件终止
  
• 垂直线的高度显示累积概率的变化
  
• 带有刻度线的经过删失的观察结果会减少间隔之间的累积生存期。
  

Kaplan-Meier图

建立在上ggplot2，并可用于创建Kaplan-Meier图。

• 默认图  带相关置信带（阴影区域）的阶跃函数（实线）。
  
• 默认情况下，显示了被检查患者的刻度线，在此示例中，该刻度线本身有些模糊，可以使用选项将其取消 censor = FALSE

估计xx年生存

生存分析中经常需要关注的一个数量是生存超过一定数量（xx）年的概率。

例如，要估算生存到11年的可能性

## Call: survfit(formula = Surv(time, status) ~ 1, data = lung)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##   365     65     121    0.409  0.0358        0.345        0.486

我们发现本研究中11年生存的机率是41％。

同时显示95％置信区间的相关上下限。

xx年生存率和生存曲线

11年存活率概率为在y轴上的点对应于11一年x轴的生存曲线。

【视频】R语言生存分析原理与晚期肺癌患者分析案例|数据分享（下）：https://developer.aliyun.com/article/1495623