生存分析

#part1 生存分析

#载入所需的R包。

library(survival) library(lubridate) library(ggsurvfit) library(gtsummary) library(tidycmprsk) library(condSURV) library(dplyr)

#载入survival包中自带的lung数据集。

lung

#为了标准统一，可以将数据集status中的1改为0，2改为1，此时，status：0=生存，1=死亡

#sex:1=男性，2=女性。

lung <- lung %>% mutate( status = recode(status, 1 = 0, 2 = 1) ) #浅看一下time, status, sex这几列。 head(lung[, c("time", "status", "sex")])

#数据通常包含开始和结束日期，而不是预先计算的生存时间。 #我们确保这些数据被格式化为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

#我们看到这两个都是字符变量，但我们需要将它们格式化为日期。 #我们将使用{lubridate}包处理日期。 #在本例中，我们需要使用ymd()函数来更改格式，日期目前是字符格式，年在前面，然后是月，最后是日。

date_ex <- date_ex %>% mutate( sx_date = ymd(sx_date), last_fup_date = ymd(last_fup_date) ) date_ex

#现在我们看到这两个日期被格式化为日期而不是字符。 #使用?ymd可以更多了解这个函数。 #现在日期已经格式化，我们需要计算开始日期和结束日期之间的差值，通常是月或年。 #使用{lubridate}包，操作符%——%指定一个时间间隔，然后使用as.duration()转换为流逝的秒数，最后通过除以dyears(1)转换为年，dyears(1)给出一年中的秒数。

date_ex <- date_ex %>% mutate( os_yrs = as.duration(sx_date %--% last_fup_date) / dyears(1) ) date_ex

#看一下前十个人的生存时间（天）。 Surv(lung status)[1:10]

#survfit()函数使用基于公式的Kaplan-Meier方法创建生存曲线，命名为s1。

s1 <- survfit(Surv(time, status) ~ 1, data = lung)

#展示一下。

str(s1)

#使用{ggsurvfit}包来生成Kaplan-Meier图。详情请参见http://www.danieldsjoberg.com/ggsurvfit/index.html。 #survfit中的time、surv用于创建生存曲线。labs函数定义x和y轴的标签。

survfit2(Surv(time, status) ~ 1, data = lung) %>% ggsurvfit() + labs( x = "Days", y = "Overall survival probability" )

#可以使用add_confidence_interval()添加置信区间。

survfit2(Surv(time, status) ~ 1, data = lung) %>% ggsurvfit() + labs( x = "Days", y = "Overall survival probability" ) + add_confidence_interval()

#可以使用add_risktable()来添加x轴下的风险数字。

survfit2(Surv(time, status) ~ 1, data = lung) %>% ggsurvfit() + labs( x = "Days", y = "Overall survival probability" ) + add_confidence_interval() + add_risktable()

#展示一下，我们发现在本研究中，1年的存活概率为41%。

summary(survfit(Surv(time, status) ~ 1, data = lung), times = 365.25) #使用{gtsummary}包中的tbl_survfit()函数生成时间生存概率估计表。 survfit(Surv(time, status) ~ 1, data = lung) %>% tbl_survfit( times = 365.25, label_header = "1-year survival (95% CI)" )

#展示中位生存时间。

survfit(Surv(time, status) ~ 1, data = lung) #使用{gtsummary}包中的tbl_survfit()函数生成中位生存时间估计表。 survfit(Surv(time, status) ~ 1, data = lung) %>% tbl_survfit( probs = 0.5, label_header = "Median survival (95% CI)" )

#log-rank检验比较组间生存时间。例如不同性别的存活时间是否存在差异？

survdiff(Surv(time, status) ~ sex, data = lung)

#我们可能想要量化单个变量的效应大小，或者在回归模型中包含多个变量，以解释多个变量的影响。 #Cox回归模型是一种半参数模型，可用于拟合具有生存结果的单变量和多变量回归模型。 #使用{survival}包中的coxph()函数为生存数据拟合回归模型。

coxph(Surv(time, status) ~ sex, data = lung)

#使用{gtsummary}包中的tbl_regression()函数获得结果表，并将取幂选项设置为TRUE以返回风险比而不是日志风险比。

coxph(Surv(time, status) ~ sex, data = lung) %>% tbl_regression(exp = TRUE)

#Cox回归模型的兴趣量是风险比(HR)。HR表示在任何特定时间点上两组之间的危险之比。

#part2 Landmark Analysis and Time Dependent Covariates

#载入所需R包

library(SemiCompRisks)

library(tidyverse)

#我们将使用{semi omprisks}包中的BMT数据集作为示例数据集。数据包括137例骨髓移植患者。感兴趣的变数包括:信息包括：T1:生存时间(以天为单位)；Delta1：1 =死亡,0 =活着； #TA:急性移植物抗宿主病的时间(天);deltaA：1=发生急性移植物抗宿主病，0=从未发生急性移植物抗宿主病。 #急性移植物抗宿主病发生在移植后90天内。 #在BMT数据感兴趣的是急性移植物抗宿主病和存活之间的关联。 #但是aGVHD是在移植后进行评估的，这是我们的基线，也就是后续随访的开始时间。

data(BMT, package = "SemiCompRisks")

#看一下前几行

head(BMT[, c("T1", "delta1", "TA", "deltaA")])

#筛出T1>=90的患者，生成新的数据lm_dat。

lm_dat <- BMT %>% filter(T1 >= 90)

#生存时间变为T1 - 90，定义为lm_T1

lm_dat <- lm_dat %>% mutate( lm_T1 = T1 - 90 )

#根据是否发生抗宿主病分组进行生存分析。

survfit2(Surv(lm_T1, delta1) ~ deltaA, data = lm_dat) %>% ggsurvfit() + labs( x = "Days from 90-day landmark", y = "Overall survival probability" ) + add_risktable()

#使用{gtsummary}包中的tbl_regression()函数查看结果。

coxph( Surv(T1, delta1) ~ deltaA, subset = T1 >= 90, data = BMT ) %>% tbl_regression(exp = TRUE)

#Part 3: Competing Risks

#我们将使用{MASS}包中的黑色素瘤数据来说明这些概念。 #status：1=死于黑色素瘤，2=活下来，3=死于其他原因；性别：1 =男性，0 =女性；ulcer：1 = 存在, 0 = 不存在。

library(MASS)

data(Melanoma, package = "MASS")

#重新定义数据集中某些变量。

Melanoma <- Melanoma %>% mutate( status = as.factor(recode(status, 2 = 0, 1 = 1, 3 = 2)) )

#看下前六行。

head(Melanoma)

#使用{tidycmprsk}包中的cuminc函数估计竞争风险背景下的累积发生率。默认情况下，这要求状态是一个因子变量，被删除的患者编码为0。

library(tidycmprsk)

cuminc(Surv(time, status) ~ 1, data = Melanoma)

#我们可以使用{ggsurvfit}包中的ggcuminc()函数来绘制累积发生率。默认情况下，它只绘制第一个事件类型。下面这张图显示的是黑色素瘤的累积死亡率。

library(ggsurvfit)

cuminc(Surv(time, status) ~ 1, data = Melanoma) %>% ggcuminc() + labs( x = "Days" ) + add_confidence_interval() + add_risktable()

#如果我们想包含这两种事件类型，请在ggcuminc(outcome =)参数中指定结果。

cuminc(Surv(time, status) ~ 1, data = Melanoma) %>% ggcuminc(outcome = c("1", "2")) + ylim(c(0, 1)) + labs( x = "Days" )

#根据溃疡，有无溃疡。我们可以按组估计不同时间的累积发病率。

cuminc(Surv(time, status) ~ ulcer, data = Melanoma) %>% tbl_cuminc( times = 1826.25, label_header = "{time/365.25}-year cuminc") %>% add_p()

#分组分析画图

cuminc(Surv(time, status) ~ ulcer, data = Melanoma) %>% ggcuminc() + labs( x = "Days" ) + add_confidence_interval() + add_risktable()

#假设我们有兴趣研究年龄和性别对黑素瘤死亡的影响，而其他原因的死亡是一个竞争事件。

#Fine-Gray回归

crr(Surv(time, status) ~ sex + age, data = Melanoma)

#我们可以使用{gtsummary}包中的tbl_regression()函数生成格式化结果表，并使用选项exp = TRUE来获得风险比估算。

crr(Surv(time, status) ~ sex + age, data = Melanoma) %>% tbl_regression(exp = TRUE)

#COX回归

coxph( Surv(time, ifelse(status == 1, 1, 0)) ~ sex + age, data = Melanoma ) %>% tbl_regression(exp = TRUE)