
生存分析(也称为工程中的可靠性分析)的目标是在协变量和事件时间之间建立联系。生存分析的名称源于临床研究,其中预测死亡时间,即生存,通常是主要目标。
生存分析是一种回归问题(人们想要预测一个连续值),但有一个转折点。它与传统回归的不同之处在于,在生存分析中,结果变量既有一个事件,也有一个与之相关的时间值,部分训练数据只能被部分观察——它们是被删失的。本文用R语言生存分析晚期肺癌患者数据
普通最小二乘回归方法不足,因为事件发生的时间通常不是正态分布的,并且模型无法处理删失,但这在生存数据中很常见。

在某些情况下,可能无法观察到事件时间:这通常称为 右删失。在以死亡为事件的临床试验中,当发生以下情况之一时,就会发生这种情况。1。当一定数量的参与者死亡时,研究结束。2。参与者退出研究。3。 研究达到预定的结束时间,并且一些参与者存活到结束。在每种情况下,幸存的参与者离开研究后,我们都不知道他们会发生什么。然后我们有一个问题:
当对于某些个体,我们只观察到他们的事件时间的下限时,我们如何对经验分布进行建模或进行非负回归?
上图说明了右删失。对于参与者 1,我们看到他们何时死亡。参与者 2 退出了,我们知道他们一直活到那时,但不知道后来发生了什么。对于参与者 3,我们知道他们活到了预定的研究结束,但又不知道之后发生了什么。

生存分析中的两个关键工具是生存函数和风险函数。
生存函数:它是一个函数,用于给出我们有兴趣知道的任何对象是否会在任何指定时间之后存活的概率。在数学上它可以由以下公式表示

其中 S(t) 是一个生存函数,其中 T 是一个连续随机变量,是一个事件的时间。F(t) 是区间[0,∞) 上的累积分布函数。
我们也可以用风险函数来写生存函数。假设事件尚未发生 ,风险率λ(t) 是事件在时间t发生的瞬时概率的主要值。


那么关键问题是如何估计风险和/或生存函数。
Kaplan Meier的非参数估计
在非参数生存分析中,我们要估计生存函数没有协变量,并且有删失。如果我们没有删失,我们可以从经验 CDF 开始. 这个等式简洁地表示:
有多少人随着时间的推移而死亡? 那么生存函数就是:还有多少人还活着?
但是,我们无法回答一些人被时间t删失时提出的这个问题.
虽然我们不一定知道有多少人在任意时间t幸存下来,我们知道研究中有多少人仍然处于风险之中。我们可以使用它来代替。将学习时间划分区间, 其中每个ti是参与者的事件时间或删失时间。假设参与者只能在观察到的事件时间失效。假设没有人在同一时间死去(没有关系),我们可以查看每次有人死去的时间。我们说在那个特定时间死亡的概率是,并说在任何其他时间死亡的概率是0.
在温和的假设下,包括参与者具有独立且相同分布的事件时间,并且删失和事件时间是独立的,这给出了一个一致的估计量。上图给出了一个简单案例的 Kaplan Meier 估计示例。

生存分析用于各种领域
例如:
在癌症研究中,典型的研究问题如下:

本演示文稿将介绍生存分析 ,参考:
Clark, T., Bradburn, M., Love, S., & Altman, D. (2003). Survival analysis part I: Basic concepts and first analyses. 232-238. ISSN 0007-0920.
我们今天将使用的一些软件包包括:
lubridatelibrary(survival)
事件时间数据由不同的开始时间和结束时间组成。
癌症的例子
事件发生时间数据在许多领域都很常见,包括但不限于
由于生存分析在许多其他领域很常见,因此也有其他名称
数据包含来自北中部癌症治疗组的晚期肺癌患者。今天我们将用来演示方法的一些变量包括
某个主题可能由于以下原因而被删失:
具体来说,这些是删失的示例。
对于主题ii:
lung数据中提供了观察时间和事件指示

数据通常带有开始日期和结束日期,而不是预先计算的生存时间。第一步是确保将这些格式设置为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我们看到它们都是字符变量,通常都是这种情况,但是我们需要将它们格式化为日期。
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-31R格式必须包含分隔符和符号。例如,如果您的日期格式为m / d / Y,则需要format = "%m/%d/%Y"我们还可以使用该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-31R选项不同,不需要指定分隔符现在日期已格式化,我们需要以某些单位(通常是几个月或几年)计算开始时间和结束时间之间的差。在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:
在lung数据中,我们有:
受试者可以存活超过指定时间的概率
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):累积分布函数
理论上,生存函数是平滑的;在实践中,我们以离散的时间尺度观察事件。
Kaplan-Meier方法是估计生存时间和概率的最常用方法。这是一种非参数方法,可产生阶跃函数,每次事件发生时,阶跃下降。
+如果主题是经过删失的,则后面跟一个。让我们看一下前10个观察值:## \[1\] 306 455 1010+ 210 883 1022+ 310 361 218 166survfit函数根据公式创建生存曲线。让我们为整个同类群组生成总体生存曲线,将其分配给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图。
plot(survfit(Surv(time, status) ~ 1, data = lung),

R中的默认图显示了具有相关置信区间(虚线)的阶跃函数(实线)建立在上ggplot2,并可用于创建Kaplan-Meier图。

censor = FALSE生存分析中经常需要关注的一个数量是生存超过一定数量(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%置信区间的相关上下限。
11年存活率概率为在y轴上的点对应于11一年x轴的生存曲线。

如果 使用“天真”的估计会怎样?
228名患者中的121名到1年时死亡,因此:

-当 忽略42名患者在1年之前受到检查的事实时, 会错误估计1个1个年生存率。

生存分析中经常需要关注的另一个数量是平均生存时间,我们使用中位数对其进行量化。预计生存时间不会呈正态分布,因此平均值不是适当的总结。
## Call: survfit(formula = Surv(time, status) ~ 1, data = lung)
##
## n events median 0.95LCL 0.95UCL
## 228 165 310 285 363我们看到中位生存时间为310天。还会显示95%置信区间的上限和下限。
中位生存时间是生存概率为0.50

总结165例死亡患者的中位生存时间
## median_surv
## 1 226lung蓝色显示,以进行比较
我们使用 函数获得对数秩p值。例如,我们可以根据lung数据中的性别测试是否存在生存时间差异
## Call:
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=1 138 112 91.6 4.55 10.3
## sex=2 90 53 73.4 5.68 10.3
##
## Chisq= 10.3 on 1 degrees of freedom, p= 0.001从 结果中提取p值
1 - pchisq(sd$chisq, length(sd$n) - 1)## \[1\] 0.001311165返回格式化的p值
## \[1\] 0.001我们可能想量化单个变量的效应大小,或者将多个变量包括在回归模型中以说明多个变量的效应。
Cox回归模型是半参数模型,可用于拟合具有生存结果的单变量和多变量回归模型。

h(t)h(t):危险或事件发生的瞬时速率h0(t)h0(t):基本基准危险
该模型的一些关键假设:
_注意_:也可以使用用于生存结果的参数回归模型,但是本培训将不涉及这些模型。
我们可以使用coxph函数拟合生存数据的回归模型,该函数Surv在左侧使用一个对象,而在右侧具有用于回归公式的标准语法R。
## Call:
##
## coef exp(coef) se(coef) z p
## sex -0.5310 0.5880 0.1672 -3.176 0.00149
##
## Likelihood ratio test=10.63 on 1 df, p=0.001111
## n= 228, number of events= 165我们可以看到输出的整洁版本broom:

或使用

estimate我们的列coxph),则HR = 经验值(β)经验值(β)。
在第1部分中,我们介绍了使用对数秩检验和Cox回归来检验感兴趣的协变量与生存结果之间的关联。
示例:从治疗开始就测量总生存期,关注的是对治疗的完全反应与生存之间的关联。
Anderson, J., Cain, K., & Gelber, R. (1983). Analysis of survival by tumor response. Journal of Clinical Oncology : Official Journal of the American Society of Clinical Oncology, 1(11), 710-9.
癌症研究中可能尚未关注的其他一些可能的协变量包括:
137例骨髓移植患者的数据。变量包括:
T1 死亡时间或最后一次随访时间(天)delta1死亡指标;1死0活TA 急性移植物抗宿主病的时间(以天为单位)deltaA急性移植物抗宿主病指标;1-发展为急性移植物抗宿主病,0-从未发展为急性移植物抗宿主病让我们加载数据以供整个示例使用
在BMT数据感兴趣的是急性移植物抗宿主病(aGVHD)和存活之间的关联。但是aGVHD是在移植后进行评估的,这是我们的基线,也就是后续随访的开始时间。
步骤1选择地标时间
通常,aGVHD发生在移植后的前90天内,因此我们使用90天的界标。
人们对急性移植物抗宿主病(aGVHD)与生存之间的关系感兴趣。但是aGVHD是在移植后进行评估的,这是我们的基线,也就是后续随访的开始时间。
第2步:至少跟踪到里程碑时间之前的人群的子集
这将我们的样本量从137减少到122。
人们对急性移植物抗宿主病(aGVHD)与生存之间的关系感兴趣。但是aGVHD是在移植后进行评估的,这是我们的基线,也就是后续随访的开始时间。
步骤3根据地标计算随访时间,并应用传统方法。

在Cox回归中, 可以使用中的subset选项coxph来排除那些在标志性时间内没有被随访的患者

界标分析的替代方法是合并时间相关的协变量。这可能更适合
对时间相关协变量的分析R需要建立特殊的数据集。
BMT数据中没有ID变量,这是创建特殊数据集所必需的,因此请创建一个名为的变量my_id。
将tmerge函数与event和函数一起使用tdc可创建特殊数据集。
tmerge 为每个患者的不同协变量值创建一个具有多个时间间隔的长数据集event 创建新的事件指示器,以与新创建的时间间隔一致tdc 创建与时间相关的协变量指标,以与新创建的时间间隔一致要了解其作用,让我们看一下前5名患者的数据。
## my_id T1 delta1 TA deltaA
## 1 1 2081 0 67 1
## 2 2 1602 0 1602 0
## 3 3 1496 0 1496 0
## 4 4 1462 0 70 1
## 5 5 1433 0 1433 0这些相同患者的新数据集
## my_id T1 delta1 id tstart tstop death agvhd
## 1 1 2081 0 1 0 67 0 0
## 2 1 2081 0 1 67 2081 0 1
## 3 2 1602 0 2 0 1602 0 0
## 4 3 1496 0 3 0 1496 0 0
## 5 4 1462 0 4 0 70 0 0
## 6 4 1462 0 4 70 1462 0 1
## 7 5 1433 0 5 0 1433 0 0现在,我们可以分析这个时间依赖性协照常使用Cox回归与coxph

我们发现,使用标志性分析或时间依赖性协变量,急性移植物抗宿主病与死亡无显着相关性。
通常,人们会希望使用地标分析对单个协变量进行可视化, 使用带有时间相关协变量的Cox回归进行单变量和多变量建模。
当对象在事件发生时间设置中发生多个可能的事件时
例子:
在任何给定的研究中,所有这些(或其中一些 以及其他)可能都是可能的事件。
事件时间之间未观察到的依赖性是导致需要特殊考虑的基本问题。
例如,可以想象复发的患者更有可能死亡,因此复发时间和死亡时间将不是独立事件。
存在多种潜在结果时的两种分析方法:
这些方法中的每一种都可能仅阐明数据的一个重要方面,而有可能使其他方面难以理解,因此所选的方法应取决于感兴趣的问题。
它包含变量:
time 生存时间以天为单位,可能经过删失。status 1例死于黑色素瘤,2例存活,3例因其他原因死亡。sex 1 =男性,0 =女性。age 年岁。year 操作。thickness 肿瘤厚度(毫米)。ulcer 1 =存在,0 =不存在。在竞争风险的背景下估算累积发生率。
## Estimates and Variances:
## $est
## 1000 2000 3000 4000 5000
## 1 1 0.12745714 0.23013963 0.30962017 0.3387175 0.3387175
## 1 3 0.03426709 0.05045644 0.05811143 0.1059471 0.1059471
##
## $var
## 1000 2000 3000 4000 5000
## 1 1 0.0005481186 0.0009001172 0.0013789328 0.001690760 0.001690760
## 1 3 0.0001628354 0.0002451319 0.0002998642 0.001040155 0.001040155生成 默认值的基本图。
plot(ci_fit)

比较组之间的累积发生率
用于组间测试。
例如,Melanoma根据ulcer溃疡的存在与否比较结果。测试结果可以在中找到Tests。
ci_ulcer\[\["Tests"\]\]## stat pv df
## 1 26.120719 3.207240e-07 1
## 3 0.158662 6.903913e-01 1
_请注意,_我个人发现该ggcompetingrisks功能缺少自定义功能,尤其是与相比ggsurvplot。我通常会自己做图,首先创建cuminc拟合结果的整洁数据集,然后再绘制结果。有关底层代码的详细信息,请参见此演示文稿的
通常,只有一种类型的事件会引起人们的兴趣,尽管我们仍要考虑竞争事件。在那种情况下,感兴趣的事件可以单独绘制。同样,我首先通过创建cuminc拟合结果的整洁数据集,然后绘制结果来手动执行此操作。有关底层代码的详细信息,请参见此演示文稿的源代码。

您可能想将风险表的数量添加到累积发生率图中,而据我所知,没有简单的方法可以做到这一点。请参阅此演示文稿的源代码中的一个示例

两种方法:
假设我们有兴趣研究年龄和性别对黑色素瘤死亡的影响,而其他原因的死亡则是竞争事件。
crr 需要指定协变量作为矩阵failcode选项请求其他事件的结果,默认情况下会返回failcode = 1shr_fit## convergence: TRUE
## coefficients:
## sex age
## 0.58840 0.01259
## standard errors:
## \[1\] 0.271800 0.009301
## two-sided p-values:
## sex age
## 0.03 0.18在上一个示例中,sex和和age均被编码为数字变量。如果存在字符变量,则必须使用model.matrix
或当前crr不支持的输出。

删失所有没有引起关注的对象,在这种情况下是由于黑色素瘤死亡,并且照常使用coxph。因此,现在对因其他原因死亡的患者进行针对特定原因的风险评估方法以应对竞争风险。


可能会出现很多零碎的东西 :
Cox比例风险回归模型的一个假设是,在整个随访过程中,风险在每个时间点都是成比例的。我们如何检查数据是否符合此假设?
使用cox.zph生存包中的功能。结果有两点:
print(cz)## rho chisq p
## sex 0.1236 2.452 0.117
## age -0.0275 0.129 0.719
## GLOBAL NA 2.651 0.266
``````
plot(cz)

有时可能想根据连续变量来可视化生存估计。求 生存数据的分位数。默认分位数是p = 0.5中位生存期。


有时,在已经存活了一段时间的患者中产生存活率估计值很有意义。

Zabor, E., Gonen, M., Chapman, P., & Panageas, K. (2013). Dynamic prognostication using conditional survival estimates. Cancer, 119(20), 3589-3592.
让我们将生存期定为6个月
map_df(
prob_times,
~conditional\_surv\_est(
basekm = fit1,
t1 = 182.625,
t2 = .x)
) %>%
mutate(months = round(prob_times / 30.4)) %>%
select(months, everything()) %>%
kable()
我们还可以根据不同的生存时间长度可视化条件生存数据。

所得出的曲线在我们每次进行条件调整时都有一条生存曲线。在这种情况下,第一条线是总体生存曲线,因为它是根据时间0进行调节的。
本文摘选《R语言中的生存分析Survival analysis晚期肺癌患者》