
本演示文稿将介绍生存分析 ,参考:
Clark, T., Bradburn, M., Love, S., & Altman, D. (2003). Survival analysis part I: Basic concepts and first analyses. 232-238. ISSN 0007-0920.
我们今天将使用的一些软件包包括:
lubridatesurvivalsurvminerlibrary(survival)library(survminer)library(lubridate)事件时间数据由不同的开始时间和结束时间组成。
癌症的例子
事件发生时间数据在许多领域都很常见,包括但不限于
由于生存分析在许多其他领域很常见,因此也有其他名称
数据包含来自北中部癌症治疗组的晚期肺癌患者。今天我们将用来演示方法的一些变量包括

RICH JT,NEELY JG,PANIELLO RC,VOELKER CCJ,NUSSENBAUM B,WANG EW。理解KAPLAN-MEIER曲线的实用指南。耳鼻咽喉头颈外科:美国耳鼻咽喉头颈外科学会官方杂志。2010; 143(3):331-336。doi:10.1016 / j.otohns.2010.05.007。
某个主题可能由于以下原因而被审查:
具体来说,这些是审查的示例。

在此示例中,我们将如何计算10年无事件的比例?
受试者2、3、5、6、8、9和10 在10年时都是无事件的。受试者4和7 在10年之前发生了该事件。主题1 在10年之前已被审查,因此我们不知道他们是否在10年之前有此事件-我们如何将该主题纳入我们的估计中?

对于主题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), xlab = "Days", ylab = "Overall survival probability")
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.266plot(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进行调节的。
参考文献
1.R语言绘制生存曲线估计|生存分析|如何R作生存曲线图
2.R语言生存分析可视化分析
3.R语言如何在生存分析与Cox回归中计算IDI,NRI指标
4.r语言中使用Bioconductor 分析芯片数据
5.R语言生存分析数据分析可视化案例
6.r语言ggplot2误差棒图快速指南
7.R 语言绘制功能富集泡泡图
8.R语言如何找到患者数据中具有差异的指标?(PLS—DA分析)
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。