生存分析R语言绘图——ggsuvplot介绍及实例

ggsurvplot(

fit, #生存分析结果

data = NULL, #a dataset used to fit survival curves

fun = NULL, # 定义生存曲线转换的任意函数。 经常使用的转换可以用字符参数指定:“event”绘制累积事件(f(y) = 1-y)，“cumhaz”绘制累积风险函数(f(y) = -log(y))，“pct”以百分比表示生存概率。

color = NULL, #曲线颜色

palette = NULL, #颜色调色板，可选调色板有 "grey","npg","aaas","lancet","jco", "ucscgb","uchicago","simpsons"和"rickandmorty".

linetype = 1, #线条形状，可以用数值型向量1，2表示，也可以用字符串向量c("solid", "dashed").

conf.int = FALSE, #是否画出置信区间

pval = FALSE, #是否显示P值

pval.method = FALSE, #是否添加计算P值得方法得文本，前提是pval = TRUE

test.for.trend = FALSE, #默认是F，如果TURE，返回trend Pvalues检验。 趋势检验旨在检测生存曲线的有序差异。 也就是说，至少对一个群体来说。 只有组数为>2时，才能进行趋势测试。

surv.median.line = "none", #画一条水平或者垂直得生存中位值线，允许的值有c("none", "hv", "h", "v"). v: 垂直vertical, h:水平horizontal.

risk.table = FALSE, #是否显示风险table。其他值有absolute" or "percentage",显示绝对数值/百分比；参数"abs_pct" ,百分比以及绝对数值都显示

cumevents = FALSE, # logical value specifying whether to show or not the table of the cumulative number of events.

cumcensor = FALSE, #logical value specifying whether to show or not the table of the cumulative number of censoring.

tables.height = 0.25, #设置table得高度，取值范围0-1

group.by = NULL, #包含分组变量名称得字符串向量。长度<=2

facet.by = NULL, #一个字符向量，包含将生存曲线分成多个面板的分组变量的名称。

add.all = FALSE, #一个逻辑值。 如果为TRUE，则在主图中添加合并患者(null model)的生存曲线。

combine = FALSE, # a logical value. If TRUE, combine a list survfit objects on the same plot.

ggtheme = theme_survminer(), #主题名称

tables.theme = ggtheme, #主题名称，默认是theme_survminer.

... #后面描述的参数和其他参数将被传递给ggplot2 geom_*()函数，如linetype, size, ii)或ggpar()函数来定制图形。 看到的细节部分

)

下图显示内置数据集colon，病人rx处理分为三组（下图第三列），对照组： Obs ，处理组一： Levamisole (Lev) ，处理组二： Levamisole + 5-fluorouracil (Lev+5FU)

# loads dplyr

library(dplyr)

# core survival analysisfunctions

library(survival)

# recommended forvisualizing survival curves

library(survminer)

#加载内置colon数据集

data(colon)

#list directory contents

ls(colon)

得到如下图：

#创建生存对象

fit <- survfit(Surv(time, status)~rx, data=colon)

#level()是为了看分组水平情况，以确定对照组，一般level（）之后第一个为对照组

levels(colon$rx)

ggsurvplot(fit,risk.table=TRUE,#生存统计统计表

                  conf.int=TRUE,#添加置信区间带

                  palette = c("skyblue","green","red"),#颜色设置

                  pval=TRUE,#log-rank检验

                   pval.method=TRUE)#添加检验text

至于是treatment中的哪一组与Obs相比，显著性，差异性更大，需要查看 Lev 和 obs 对比的p值及HR，以及 （Lev+5FU） 和 Obs 对比的p值及HR，评价分组的治疗效果

#Cox Regression,评价rx分组后治疗效果

fit1<-coxph(Surv(time, status)~rx, data=colon)

fit1

我们经常用随机森林等机器学习又或者是其他数据挖掘的方法寻找某些疾病的biomarker或者候选基因。但是来自临床的数据包括了生存事件等信息，数据的内容有所不同，所以需要一些和之前不太一样分析方法，其中常见的就是通过制作生存曲线图获取结论。

生存曲线可以帮助我们回答许多问题：参与者生存5年的概率是多少？两组之间的生存率是否存在差异（例如，在临床试验中分配给新药还是标准药的两组之间）？某些行为或临床特征如何影响参与者的生存机会？

通常，在这类分析中，我们会关注特定事件（如死亡或疾病复发）的事件，并比较两组或更多组患者发生这些特定事件的事件。

可以看到上图显示了经常玩棋类游戏的老年人和很少玩这类游戏的老年人之间的痴呆风险Kaplan-Meier曲线。纵轴为非痴呆老人的比例，横轴为跟踪的年数，从图中可以看到经常玩棋类游戏的老年人患痴呆的风险较低。

在制作生长曲线之前，我们需要首先了解几个相关的术语

参考： R语言-Survival analysis（生存分析）

Event（事件）： 指在随访过程中发生的某个结果，如癌症研究中，可能为复发（Relapse）、恶化（Progression）、死亡（Death）

Survival time（生存时间）： 指某个事件开始到终止的时间，在癌症研究中经常用到的几个指标：

Overall survival（OS）：

指从开始到任意原因死亡的时间，一般常见的5年生存率、10年生存率都是基于OS计算的

Progression-free Survival（PFS，无进展生存期）：

指从开始到肿瘤发生任意进展或者死亡的时间，可用于评估治疗方法的临床效益

Time to Progress（TTP，疾病进展时间）：

从开始到肿瘤发生任意进展或者进展前死亡的时间，与PFS相比仅包括肿瘤的恶化，而不包括死亡。

Disease-free Survival（DFS，无病生存期）：

指从开始到肿瘤复发或任何原因死亡的时间，常用于根治性手术治疗或放疗后的辅助治疗的评估

Event Free Survival（EFS，无事件生存期）：

指从开始到发生包括肿瘤进展、死亡、治疗方案的改变等各种事件的时间

Censoring（删失）： 一般指不是由于死亡造成的数据丢失，可能是由于失访、非正常原因推出、时间终止而事件未发生等，一般在展示时用“+”表示

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

1、参数法：已知生存时间的分布模型，根据数据估计模型参数，最后以分布模型计算生存率

2、半参数法：不需要知道生存时间的分布，但是仍通过模型来评估影响生存率的因素，常见方法如 Cox回归模型

3、非参数法：不需要知道生存时间的分布，根据样本统计量估计生存率，常见方法如 Kaplan-Meier方法、寿命法

具体地，我们通过同样一个例子介绍常用的Kaplan-Meier方法和寿命法的异同。

例子：一项探究死亡时间的前瞻性队列研究，研究涉及20位65岁以上的参与者，招募时间为5年，整个研究进行长达24年的随访直至死亡、研究结束或退出研究（失访）。因此，如果参与者是在研究开始后加入的，他们的最长随访时间应该少于24年。具体数据如下，其中有6位参与者死亡，3位接受了完整的随访（24年），其余11位由于在研究开始后加入或失访而少于24年随访:

寿命法

寿命法经常用于保险行业中估计预期寿命并设置保费。不过，我们只关注生物领域的使用，我们称为随访生命表，该表记录了参与者在队列研究或临床试验中在预定的随访期内的经历，直到目标事件发生或研究结束为止。

要构建生命表，我们要将随访时间分割成间距相等的几组，上述例子中我们随访的最长时间为24年，所以我们考虑5年一个间隔（0-4，5-9，10-14，15-19和20-24年）。然后统计每个时间间隔开始时活着的参与者人数，和该期间死亡人数和每个时间间隔中删失的人数。

然后，我们来定义几个参数：

N t =在时间间隔t内没有发生目标事件的但处于风险中的人数（如本研究中目标事件为死亡，而参与者都处于可能死亡的风险之中）

D t =在时间间隔t内死亡的人数

C t =在时间间隔t内删失的人数

N t * =在时间间隔t内有风险的参与者的平均数（计算公式为：N t * =N t -C t /2）

q t =时间间隔t内死亡比例，q t =D t /N t *

p t =时间间隔t内生存比例，p t =1-q t

S t ，累计生存概率，S 0 =1,S t+1 =p t+1 *S t

因此，对于第一个间隔0-4年和第二个5-9年的间隔，可以计算出如下数据：

所以完整的随访寿命表为：

Kaplan-Meier

Edward Kaplan和Paul Meier于1958年在《American Statistical Association》共同发表了Kaplan-Meier非参数估计方法，让我们能够估计生存函数。

从寿命表的方法可以看出生存概率会根据不同的间隔改变，尤其是对于小样本而言这种改变可能会很剧烈。Kaplan-Meier通过每次时间发生时重新估计生存概率来解决该问题。

Kaplan-Meier是基于这样的假设进行的：删失与事件发生的可能性无关，且在研究早期和后期被招募的参与者生存率是可比的。这些前提很重要，比如在不同组比较时要保证删失的可能性一致。

Kaplan-Meier与寿命法的计算方式类似，主要区别是时间间隔，寿命法中我们选择的时间间隔相等，而在Kaplan-Meier的方法中我们使用观察到的事件时间和删失时间。

上述的内容原版，以及关于进一步的检验和Cox模型的内容可以阅读Boston大学的教材 Boston Univeristy Suvival Analysis 。在这里暂时就不再解释啦。

今天我们要用到以下几个R包：survival，survminer和dplyr

使用KM方法，通过 ggsurvplot 作图，该函数作图需要两部分数据，具体见下：

1）需要什么格式的数据

我们使用的数据集为ovarian，来自survival包。该数据集来源于文章《Different Chemotherapeutic Sensitivities and Host Factors Affecting Prognosis in Advanced Ovarian Carcinoma vs. Minimal Residual Disease》，主要研究化疗敏感性和宿主因素对晚期卵巢癌和微小残留病变的预后影响，具体含有以下几个指标：

futime: survival or censoring time 生存时间

fustat: censoring status 确定参与者生存时间是否发生缺失

age:in years

resid.ds: residual disease present (1=no,2=yes) 评估肿瘤的消退情况

rx: treatment group 接受两种治疗方案中的一种

ecog.ps:ECOG performance status (1 is better, see reference)依据ECOG评估的患者表现

为了更直观的获取信息，我们根据说明修改一下部分指标的标记方式：

然后我们来看一下年龄的分布 hist(ovarian$age)

然后我们进行生存曲线的分析，使用futime和fustat两列，首先根据是否发生删失对数据进行处理。

可以看到发生删失的都带上了加号。

然后拟合Kaplan-Meier曲线：

2)如何作图

然后使用 ggsurvplot 功能进行绘图，如果选择 pval=TRUE 会显示两组差异检验结果的pvalue。

如果想要研究与resid.ds的关系：

往期R数据可视化分享

R数据可视化13:瀑布图/突变图谱

R数据可视化12: 曼哈顿图

R数据可视化11: 相关性图

R数据可视化10: 蜜蜂图 Beeswarm

R数据可视化9: 棒棒糖图 Lollipop Chart

R数据可视化8: 金字塔图和偏差图

R数据可视化7: 气泡图 Bubble Plot

R数据可视化6: 面积图 Area Chart

R数据可视化5: 热图 Heatmap

R数据可视化4: PCA和PCoA图

R数据可视化3: 直方/条形图

R数据可视化2: 箱形图 Boxplot

R数据可视化1: 火山图

---