之前写过生存分析的数学相关基础知识，这次直接使用R语言进行生存分析的实战演练。

1. 生存分析

install.packages(c("survival", "survminer"))
# Load the packages
library("survival")
library("survminer")

导入示例需要的数据

# Example data sets
data("lung")
head(lung)
dim(lung)
colnames(lung)

看一看这个数据里面都有什么

inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
1    3  306      2  74   1       1       90       100     1175      NA
2    3  455      2  68   1       0       90        90     1225      15
3    3 1010      1  56   1       0       90        90       NA      15
4    5  210      2  57   1       1       90        60     1150      11
5    1  883      2  60   1       0      100        90       NA       0
6   12 1022      1  74   1       1       50        80      513       0

有这几项数据

> colnames(lung)
 [1] "inst"      "time"      "status"   
 [4] "age"       "sex"       "ph.ecog"  
 [7] "ph.karno"  "pat.karno" "meal.cal" 
[10] "wt.loss"  

这几项数据的含义如下：

inst: Institution code 机构代码
time: Survival time in days 生存时间
status: censoring status 1=censored, 2=dead 生存状态
age: Age in years 年龄
sex: Male=1 Female=2 性别
ph.ecog: ECOG performance score (0=good 5=dead)
ph.karno: Karnofsky performance score (bad=0-good=100) rated by physician
pat.karno: Karnofsky performance score as rated by patient ，Karnofsky表现评分，按患者评分
meal.cal: Calories consumed at meals 餐时消耗的卡路里
wt.loss: Weight loss in last six months 最后6个阅体重损失量

上面的数据，其中做单纯的生存分析，重要的就是生存时间，生存状态，基于这两项，有分组相关的其他指标，如果按照年龄划分，第三个变量就是年龄，如果按照其他分组指标就按照其他指标衡量。如前面所说，KM生存分析是一个参数检验模型，如果有多个变量，那么可以采用Cox回归分析。

2. 计算生存曲线

这个例子中利用性别计算生存率
可以使用survival package中的survifit()函数来计算KM生存计算。主要包括两个方面的内容：

• 使用surv()函数形成的生存对象
• 包含变量的数据
  针对以上的数据，做以下计算

fit <- survfit(Surv(time, status) ~ sex, data = lung)
print(fit)

结果如下

Call: survfit(formula = Surv(time, status) ~ sex, data = lung)
        n events median 0.95LCL 0.95UCL
sex=1 138    112    270     212     310
sex=2  90     53    426     348     550

默认状态下，使用print()函数只是对生存曲线计算的一个简要总结，只是展示观察数量，事件发生数量，中位生存期和中位数的置信限，如果要得到更为详细的生存曲线数据，可以使用下面的函数:

# Summary of survival curves
summary(fit)
# Access to the sort summary table
summary(fit)$table

3. 查看函数survfit()返回的值

函数survfit（）返回变量列表，包括以下数据：

n: total number of subjects in each curve. 每条曲线中的对象总数
time: the time points on the curve. 曲线上的时间点
n.risk: the number of subjects at risk at time t。 在t时刻有风险的受试者人数
n.event: the number of events that occurred at time t. 在t时刻发生的事件数
n.censor: the number of censored subjects, who exit the risk set, without an event, at time t. 在t时刻无事件退出风险集的受检主体数量
lower,upper: lower and upper confidence limits for the curve, respectively. 曲线的置信区间
strata: indicates stratification of curve estimation. If strata is not NULL, there are multiple curves in the result. The levels of strata (a factor) are the labels for the curves. 表示曲线估计的分层。

可以按照如下的方式，查看这些数据

d <- data.frame(time = fit$time,
                  n.risk = fit$n.risk,
                  n.event = fit$n.event,
                  n.censor = fit$n.censor,
                  surv = fit$surv,
                  upper = fit$upper,
                  lower = fit$lower
                  )
head(d)

结果如下

 time n.risk n.event n.censor      surv     upper     lower
1   11    138       3        0 0.9782609 1.0000000 0.9542301
2   12    135       1        0 0.9710145 0.9994124 0.9434235
3   13    134       2        0 0.9565217 0.9911586 0.9230952
4   15    132       1        0 0.9492754 0.9866017 0.9133612
5   26    131       1        0 0.9420290 0.9818365 0.9038355
6   30    130       1        0 0.9347826 0.9768989 0.8944820

4. 可视化生存数据

可以使用survminer package 中的ggsurvplot()函数对生存数据分析的可视化。这个函数能够画出分成两组的生存曲线。

# Change color, linetype by strata, risk.table color by strata
ggsurvplot(fit,
          pval = TRUE, conf.int = TRUE,
          risk.table = TRUE, # Add risk table
          risk.table.col = "strata", # Change risk table color by groups
          linetype = "strata", # Change line type by groups
          surv.median.line = "hv", # Specify median survival
          ggtheme = theme_bw(), # Change ggplot2 theme
          palette = c("#E7B800", "#2E9FDF"))

得到的生存曲线是这个样子，如果没有定义相应的颜色，默认的颜色是红蓝色，大部分文章都用这种颜色。

上面的代码中颜色、线条都是按照分组进行改变的。
当然上面的图形依然是可以进行改变的。

ggsurvplot(
   fit,                     # survfit object with calculated statistics.
   pval = TRUE,             # show p-value of log-rank test.
   conf.int = TRUE,         # show confidence intervals for 
                            # point estimaes of survival curves.
   conf.int.style = "step",  # customize style of confidence intervals
   xlab = "Time in days",   # customize X axis label.
   break.time.by = 200,     # break X axis in time intervals by 200.
   ggtheme = theme_light(), # customize plot and risk table with a theme.
   risk.table = "abs_pct",  # absolute number and percentage at risk.
  risk.table.y.text.col = T,# colour risk table text annotations.
  risk.table.y.text = FALSE,# show bars instead of names in text annotations
                            # in legend of risk table.
  ncensor.plot = TRUE,      # plot the number of censored subjects at time t
  surv.median.line = "hv",  # add the median survival pointer.
  legend.labs = 
    c("Male", "Female"),    # change legend labels.
  palette = 
    c("#E7B800", "#2E9FDF") # custom color palettes.
)

经过改变之后，得到

上面的图形都是什么含义：
横轴表示时间，纵轴表示生存的可能性或生存的人口比例。曲线中的垂直下降表示事件的发生。

• 在0时刻，生存概率为1
• 在250天，对于两种性别生存概率分别是0.55和0.75 。
• 对于性别1，中位生存时间大约为270天，性别2的中位生存时间为426天，意味着，相比于性别1，性别2有着较好的生存结果。

每一组的中位生存时间可以通过以下方式获得

summary(fit)$table

结果如下

   records n.max n.start events   *rmean *se(rmean) median 0.95LCL 0.95UCL
sex=1     138   138     138    112 325.0663   22.59845    270     212     310
sex=2      90    90      90     53 458.2757   33.78530    426     348     550

同样可以画出累计风险的图形

ggsurvplot(fit,
           conf.int = TRUE,
           risk.table.col = "strata", # Change risk table color by groups
           ggtheme = theme_bw(), # Change ggplot2 theme
           #palette = c("#E7B800", "#2E9FDF"),
           fun = "event",
           pval = T)

对于图中p值的位置，可以通过AI后期调整，也可以改变参数进行调整。参数调整的方式以后另更。

性别= 1（男性）的中位生存时间为270天，而性别= 2（女性）则为426天。
与男性相比，女性肺癌似乎具有生存优势。
但是，要评估这种差异是否具有统计学显着性，需要进行正式的统计检验，这将在下面中讨论。

• 注意，置信极限在曲线的尾部很宽，很难进行有意义的解释。
  这可以通过以下事实来解释：在实践中，通常有一些患者在随访结束时迷失了随访或存活。
  因此，在随访结束之前在x轴上缩短图可能是更合理的。

所以下面的图形就是将横轴进行了缩减得到。

ggsurvplot(fit,
           conf.int = TRUE,
           risk.table.col = "strata", # Change risk table color by groups
           ggtheme = theme_bw(), # Change ggplot2 theme
           #palette = c("#E7B800", "#2E9FDF"),
           xlim = c(0, 600),
           pval = T,
           break.time.by = 200)

image.png

5. Kaplan-Meier生命表：生存曲线摘要

上面谈及到可以用summary()函数获得生存曲线完整的摘要。同样也可以使用survminner package中的surv_summary() 函数得到生存曲线的摘要。相比于默认的summary()函数，surv_summary()能够创建一个数据框

res.sum <- surv_summary(fit)
head(res.sum)

结果如下

> View(res.sum)
> head(res.sum)
  time n.risk n.event n.censor      surv    std.err     upper     lower strata sex
1   11    138       3        0 0.9782609 0.01268978 1.0000000 0.9542301  sex=1   1
2   12    135       1        0 0.9710145 0.01470747 0.9994124 0.9434235  sex=1   1
3   13    134       2        0 0.9565217 0.01814885 0.9911586 0.9230952  sex=1   1
4   15    132       1        0 0.9492754 0.01967768 0.9866017 0.9133612  sex=1   1
5   26    131       1        0 0.9420290 0.02111708 0.9818365 0.9038355  sex=1   1
6   30    130       1        0 0.9347826 0.02248469 0.9768989 0.8944820  sex=1   1

查看数据框结果如下

6. 使用log-rank检验比较生存曲线

log-rank检验是比较两个或者多个生存曲线最常用的方法。假设为两组之间没有统计学的差异。log-rank 检验是非参数检验，不需要对生存分布做相应的假设。
本质上，log-rank检验将在每个组中观察到的事件数与如果原假设为真（即，生存曲线相同）所期望的事件数进行比较。

survival package中的survdiff()函数可以用来计算log-rank检验中比较两个或更多生存曲线。方法如下：

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

结果为

Call:
survdiff(formula = Surv(time, status) ~ sex, data = lung)
        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.00131 

检验结果p= 0.00131 ，表明按照性别分组在生存方面有统计学差异。

这篇文章主要叙述了单因素生存分析的相关计算和分析，如前一篇文章说到如果有多元统计模型建模需求，那么就需要之谈到的cox回归模型，下一次再谈用R计算cox比例风险模型。

参考文章

Survival Analysis Basics
Cox Proportional-Hazards Model
R|生存分析 - KM曲线 ，值得拥有姓名和颜值