adobe仿宋std r常规_R生存分析

概述

• Kaplan-Meier曲线描述了生存函数。
• 截尾是一种生存分析特有的缺失数据问题。
• 比例风险假设：生存分析的主要目标是比较不同组别（例如白血病患者与正常对照组）的生存功能。
• Censoring（删失）：这经常会在临床资料中看到，生存分析中也有其对应的参数，一般指不是由死亡引起的的数据丢失，可能是失访，可能是非正常原因退出，可能是时间终止而事件未发等等，一般在展示时以‘+’号显示
  • left censored（左删失）：只知道实际生存时间小于观察到的生存时间
  • right censored（右删失）：只知道实际生存时间大于观察到的生存时间
  • interval censored（区间删失）：只知道实际生存时间在某个时间区间范围内

方法分类

生存分析的方法一般可以分为三类：

• 参数法：知道生存时间的分布模型，然后根据数据来估计模型参数，最后以分布模型来计算生存率
• 非参数法：不需要生存时间分布，根据样本统计量来估计生存率，常见方法Kaplan-Meier法（乘积极限法）、寿命法
• 半参数法：也不需要生存时间的分布，但最终是通过模型来评估影响生存率的因素，最为常见的是Cox回归模型

比例风险回归也称为Cox回归，是评估不同变量对生存率影响的最常用方法。

R包survival和survminer

生存分析的R包一般用survival包和survminer包。核心的分析函数都在survival包里，然后用survminer包进行绘图。

survival核心函数包括：

• Surv()：创建一个生存对象
• survfit()：使用公式或已构建的Cox模型拟合生存曲线
• coxph()：拟合Cox比例风险回归模型

另外有函数：

• cox.zph()：检验一个Cox回归模型的比例风险假设
• survdiff()：用log-rank/Mantel-Haenszel检验检验生存差异

Surv()创建响应变量，典型用法是使用事件时间，以及事件是否发生（即死亡与截尾）。 survfit()创建一条生存曲线，然后可以显示或绘图。 coxph()实现回归分析，并且模型以与常规线性模型（ lm()、 glm()）中相同的方式指定。

data source

library(survival)
data(lung)
str(lung)
# 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
# meal.cal: Calories consumed at meals
# wt.loss: Weight loss in last six months

Surv对象构建

函数：

Surv(time, time2, event, type=c('right', 'left', 'interval', 'counting', 'interval2', 'mstate'),origin=0)

• time: 对于右删失，即指随访的时间。若为区间数据，则为区间数据的开始时间；
• time2: 区间数据，则为区间数据的结束时间；
• event: 制定结局变量。默认0为删失，1为出现终点事件。可以自己指定终点事件。df$status == 1表示认为值为2时为TRUE，即提示终点事件发生。同理df$status == 2提示终点事件发生
• type:制定删失类型。但是type在不特别指定的情况下，一般会自动根据数据进行默认选择。所以不用单独进行指定，默认值就很好。

surv <- Surv(time = lung$time, lung$status == 2)
# lung$surv <- Surv(time = lung$time, lung$status == 2)

• surv的结果只是一列一连串的生存时间汇总。所以一般会提倡用lung$surv这种方式把结果整合到原始数据框中便于之后的分析

生存曲线与结果汇总

使用survfit()函数拟合一条生存曲线。如果要创建一条不考虑任何比较的生存曲线，只需要指定survfit()在公式里期望的截距为~1。这一点与lm()、glm()相似。

survfit(formula = lung$surv~1)
# Call: survfit(formula = s ~ 1)
# 
# n  events  median 0.95LCL 0.95UCL 
# 228     165     310     285     363

模型对象本身不会给出太多的价值信息，我们需要使用summary函数查看模型汇总结果

summary(survfit(formula = lung$surv~1, data = lung))
# 表格每一行显示了一个（多个）事件或截尾发生了，
# 在风险中的样本数（就是还没死的），以及及时的累积生存率等
# * time: the time points on the curve.
# * n.risk: the number of subjects at risk at time t
# * n.event: the number of events that occurred at time t.
# * n.censor: the number of censored subjects, who exit the risk set, without an event, at time t.
# * surv: estimate of survival probability.
# * std.err: standard error of survival.
# * upper: upper end of confidence interval
# * lower: lower end of confidence interval
# * 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.

如果加入其实分组变量，则summary返回分组条件下的结果

sfit <- survfit(formula = surv~sex, data = lung)
summary(sfit)

summary()函数中可以设定时间参数用来选定一个时间区间

summary(sfit, times=seq(0, 1000, 100))

Kaplan-Meier曲线

library(survminer)
ggsurvplot(sfit, data = lung)
# ggsurvplot() is a generic function to plot survival curves. 
# Wrapper around the ggsurvplot_xx() family functions. 
# Plot one or a list of survfit objects as generated by the survfit.formula() and surv_fit functions:
# ggsurvplot_list() 绘制多个对象
# ggsurvplot_facet() 分面到多个panels
# ggsurvplot_group_by() 一幅图中多个分组
# ggsurvplot_add_all() 总合所有的情况
# ggsurvplot_combine() 一个图中结合多个survfit对象

还可以添加曲线的置信区间，并增加long-rank检验的结果p值以及风险表格

ggsurvplot(sfit, conf.int=TRUE, pval=TRUE, risk.table=TRUE, 
           legend.labs=c("Male", "Female"), legend.title="Sex",  
           palette=c("dodgerblue2", "orchid2"), 
           title="Kaplan-Meier Curve for Lung Cancer Survival", 
           risk.table.height=.15)

分组检验

survdiff(surv~sex, data = lung, rho = 0)
# With rho = 0 this is the log-rank or Mantel-Haenszel test
# and with rho = 1 it is equivalent to the Peto & Peto modification of the Gehan-Wilcoxon test.

COX回归模型

之所以要用到COX回归模型，是因为Kaplan-Meier曲线用来对二分类变量差异的可视化比较合适，但当分类比较多时，如好、中、差、极差，效果不太理想；而且Kaplan-Meier曲线对连续型变量的风险的可视化不支持。

Cox回归模型是处理多分类变量问题的优秀方法，并且能够在对连续型变量截断后进行有效评估。其内置于survival包中，语法与lm()和glm()一致。

COX 回归模型不用于估计生存率，主要用于因素分析。

首先与KM方法中sex因素条件下的风险模型进行比较：

fit <- coxph(surv~sex, data=lung)
summary(fit)
# Call:
# coxph(formula = s ~ sex, data = lung)
# 
# n= 228, number of events= 165 
# 
# coef exp(coef) se(coef)      z Pr(>|z|)   
# sex -0.5310    0.5880   0.1672 -3.176  0.00149 **
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# exp(coef) exp(-coef) lower .95 upper .95
# sex     0.588      1.701    0.4237     0.816
# 
# Concordance= 0.579  (se = 0.021 )
# Likelihood ratio test= 10.63  on 1 df,   p=0.001
# Wald test            = 10.09  on 1 df,   p=0.001
# Score (logrank) test = 10.33  on 1 df,   p=0.001

结果中的exp(coef)列包含$e^{β_1}$。它就是风险比率——该变量对风险率的乘数效应（对于该变量每个单位增加的）。因此，对于像性别这样的分类变量，从男性（基线）到女性的结果大约减少约40％的危险。你也可以翻转coef列上的符号，并采用exp(0.531)，你可以将其解释为男性导致危险增加1.7倍，或者单位时间男性的死亡率约为女性1.7倍（女性死亡率为男性的0.588倍）。

要记住： HR=1: 没有效应 HR>1: 风险增加 * HR<1: 风险减少 （保护变量） 还能注意到，“sex”有一个对应的p值，整个模型中也有一个p值。0.00111的p值非常接近我们在Kaplan-Meier图上看到的p=0.00131的p值。这是因为KM曲线显示的是对数秩检验的p值，通过调用summary(fit)来获得Cox模型结果。你也可以使用survdiff()直接计算log-rank测试p值。

ggforest()里的图提示了HR的大小，可以告知因素是危险因素还是有益因素

cox回归模型假设检验

对coxph()的结果应该应用cox.zph()进行检验，观察global的p值，判定是否符号cox回归假设。用plot(cox.zph())来可视化每一个变量在cox模型下的表现。当数据出现明显分离时，则说明表现不是很理想。比如下图：

cox回归变量筛选

survival包中进行cox回归分析的主要函数是coxph，该函数和R中其他建模函数的用法一样，使用因变量~自变量的形式，分类变量同样是以最小编码值为参照水平。coxph函数可以使用逐步回归，采用AIC准则进行变量的筛选。下面的代码中我们首先建立空模型（仅仅包含截距项的模型），然后建立全模型（包含所有自变量的模型），进而采用三种方法进行变量选择。

其变量筛选的流程与lm()或者glm()相似，先对所有变量所回归，看p值，再用leaps包中的regsubsets()全子集回归预测，用Plot()作图，结合car包中的subsets()做全子集回归，然后作图。最后根据回归结果筛选变量，再次进行回归，对不同的回归模型进行比较，用到anova()或者chisq.test()。最后根据AIC原则纳入变量进行最后的分析。

基于cox回归模型的预测

得到cox回归模型的model对象后，可以根据model对新的数据进行预测。这些预测可以包括生存时间，生存率等等。其执行方式与机器学习部分的内容相似。主要用到predict()函数

cut()连续变量

lung数据集中age是连续型变量，首先直接用coxph()进行回归分析，在得到确证该变量为重要变量后，再进一步对其进行cut后讨论年龄的分界因素如何。比如，多少以上或者以下对生存率是如何影响的。

为什么不直接用survfit()呢？

因为直接用survfit()的话，它会把所有的年龄当作一个factor，然后给出所有factor条件下的KM曲线，那将是很多很多线的交织，但却没有任何意义。因此，将连续变量cut后，可以对其进行二分类，或者三分类，再根据对其进行的分类新变量数据采用survfit()评估风险，最后得到的KM曲线也将只基于这些分类了。

lung$age_cut <- cut(lung$age, breaks = c(0, 70, Inf), labels = c("young", "old"))
fit <- survfit(surv~age_cut, data = lung)
summary(fit)
ggsurvplot(fit, data = lung) # 不要忘记data = lung

有一点特别需要注意的，为什么年龄选择在70岁做为分界呢？这里其实涉及到分位数的概念。当我们取平均值的时候，KM曲线并未显示出差别来，但是当我们取四分位上四分之一，即70岁时，KM曲线在两个年龄层即显示出了差异。
还要记住，Cox回归分析整个分布范围内的 连续变量，其中Kaplan-Meier图上的对数秩检验可能会根据您对连续变量进行分类而发生变化。他们以一种不同的方式回答类似的问题：回归模型提出的问题是“年龄对生存有什么影响？”，而对数秩检验和KM图则问：“那些不到70岁和70岁以上的人有差异吗？“

survfit()与coxph()的作图说明

根据survminer包的描述，ggsurvplot()的作图对象是survfit()的结果，而ggforest()的作图对象是coxph()的结果。因为前者画的是生存曲线，而后者是对各变量因素风险的评估。两者可视化意义是不一样的。

如果要强行对coxph()的结果进行ggsurvplot()绘图，则需要加survfit()。不过，需要理解的是，此时coxph()的结果对survfit()来说更像只是一种表达式的存在，并没有什么意义。

所以压根儿就不要干这种事儿。 但有一种情况是特殊的，即在对新的数据进行预测时则是有意义的。不过，这个时候的 coxph()结果实际上是预测用的 model的存在， survfit()需要其中的 model去预测 新数据的KM曲线情况。一般情况下， coxph()中的变量因素会比较多。甚至会涉及到 strata()。

那么strata()到底是什么呢？即分层变量，它用于分层分析，能够控制混杂因素。若选入该变量，结果将按该变量的各水平分别输出。因为在不同的变量往往可能会有相关性，我们需要排除这种相关性的混杂因素，去纯粹地提炼出我们感兴趣的那个因素对生存率的影响。比如：

cox.test <- coxph(surv~age_cut + strata(sex), data = lung)
ggsurvplot(survfit(cox.test), data = lung)

这里就考虑了sex与age_cut的相关性。我们需要在校正了sex对生存率的影响后，弄清楚年龄是不是仍然是有意义的因素。这个时候，你在绘制图中看到两条线，分别属于两个性别，但是如果你加入pval = TRUE这个选项，你是得不到统计检验p值的。因为你无法比较的sex类别，这里它只是你希望去除的混杂因素。当然，即使你去除strata(sex)，再继续作图，并加入pval = TRUE，你还是不会得到cut年龄后的两分类的差异。因为coxph()的方法并不是用来检测同一变量内的差异的，而是用来检验该因素是否有风险的。它只会在你用ggsurvplot()调用coxph()结果时，给你一个all的总曲线，而不会自动给你分层处理，所以并没意义。

所以，本质上，你要理解coxph()和survfit()它们分别是干什么的。这是两种方法，所为不同的目的。

那么coxph()结果应该用什么来可视化呢？当然是ggforest()。ggforest()用于对coxph()结果可视化，正如其函数介绍所说。用ggforest()可视化coxph()结果，尤其适用于对包含很多个因素进行同时评价时。想象一下森林图的样子，那就很直观的显示哪些有意义，哪些没有意义了。不过，要注意的是，在进行整体因素评估时，要在formula=.-time-status，即剔除用于构建Surv()对象进的这两个因素，否则会影响结果的准确性。

coxfits <- coxph(surv~age_cut+sex+ph.ecog+ph.karno+wt.loss, data = lung)
coxfits_c <- coxph(surv~.-time-status, data = lung)
ggforest(coxfits)
ggforest(coxfits_c)

分层比较

层（Strata）：分层变量，用于分层分析，控制混杂因素。若选入变量，结果将按该变量的各水平分别输出。

survfit()函数还可以进一步进行生存率的分层比较，只需要在survfit()参数中增添层参数strata()就可以了。在survdiff()假设检验中同样如此，如下：

survdiff(surv~age_cut+strata(ph.ecog),data = lung)

多重插补

我注意到用complete.cases()判断lung数据集时，其行缺失高达67个，而总共obs才228。根据回归模型的要求，在进行回归的过程中是不允许存在NA值的。而回归模型在进行回归时，会自动删除掉带有缺失值的行。这会造成大量的obs的运算性删除。虽然有的人可能会说，如果选某个单个变量进行回归分析时，可能影响就没有那么大。但是根据指南介绍的，如果不进行插补，那么每一个回归运算由于缺失了不同行的数据，实际上运用的是不同的数据集，其结果是不理想的。所以在此建议用mice包进行多重插补后再运算。

library(survival)
library(survminer)
data(lung)
sum(!complete.cases(lung))
# [1] 67
md.pattern(lung, 5)
lung_cmplt <- mice(lung, 5) # 5种可选插补
lung_cmplt_3 <- complete(lung_cmplt, 3) # 选取其中一个完整的子集进行后续分析
lung_cmplt_3$surv <- Surv(lung_cmplt_3$time, lung_cmplt_3$status == 2)
fit <- survfit(surv~age, data = lung_cmplt)
ggsurvplot(fit, pval = TRUE)

其实我不太确定这样人为的选择其中的一个插补集进行分析是否合理。因为就常理讲，mice包在分析 lm()或者 glm()时，其分析对象就是lung_cmplt，自然也包括其中的5个可选集，最后在回归分析时，也会考量这5个可选集各自的情况，然后给出总的结论。而如果我只是选一个的话，就不存在汇总多种情况的步骤了。因此可能影响对结果的判断。或许对这几个子集进行选分析，最后再决定选哪个要比一开始主观的选择更合理一点。