生存分析如何绘制事件发生累计概率曲线图？

在生存分析中我们通常关注个体在时间t之前存在的概率(或者不发生感兴趣事件的概率)，并且绘制生存函数曲线。

生存函数(S(t)=Pr(T>t))是一个阶梯函数，说明随着时间的推移累计生存概率。但有时候研究者更加关注的是随着时间推移事件发生累计概率(顾名思义就是1-S(t))。

近期就收到一个咨询，希望可以绘制1-S(t)曲线，并且输出风险表(risk table)，在图上标注曲线的置信区间。

我们知道SPSS可以输出生存函数、1-生存函数、风险函数、对数风险函数。

SPSS具体操作如下：

1.打开SPSS软件，导入需要分析的数据，数据设置好格式后，点击分析–生存分析–Kaplan-Meier，如下图1所示：

2.如下图2所示定义好时间、状态和因子，点击右边的选项；全部勾选生存函数、1减去生存函数、风险函数、对数生存。

3.选好之后，点击继续–确定就可以得到想要的图。

其中第二张图就是累计事件发生概率曲线(1-S(t))。

但是SPSS输出图不够美观，并且不能添加置信区间和风险表。

这个时候首先想到用R的”survminer”包实现优雅地绘制生存曲线，解决这个需求最称手的当属ggsurvplot函数。

话不多说直接上代码和结果：

R软件具体操作如下：

1.导入数据(略)，加载需要的包：

library(ggplot2) # 加载ggplot2包
library(survival) # 加载survival包拟合生存对象
library(survminer) # 加载survminer包绘制生存曲线

2.定义需要用于绘制生存函数的formula，直接用ggsurvplot画图。

ggsurvplot(survfit(Surv(Os,Status)~Group, data=data1),
           pval=TRUE,
           risk.table = TRUE)
# 参数解释
Surv(Os,Status)~Group,data=data1 # 生存函数formula
pval=TRUE # 输出log-rank检验P值
pval.method=TRUE  # 可以显示检验方法
risk.table = TRUE # 输出风险表，默认FALSE

上面的代码与下面2行代码等同：

f <- survfit(Surv(Os,Status)~Group,data=data1)
ggsurvplot(f, pval=TRUE, risk.table = TRUE)

输出的图如下：

我们发现默认输出的图是生存函数图，那要怎么得到事件发生的累计概率曲线图呢？

这时候需要用到(fun = ” “)这个参数。

小编查阅了知乎、简书等分享的文章，其中不乏爱好R语言的临床研究者梳理的生存分析操作流程，有一些误以为fun="cumhaz"（累计风险图）就是事件发生的累计概率曲线图。

ggsurvplot(survfit(Surv(Os,Status)~Group,data=data1),
           fun="cumhaz",
           pval=TRUE,
           risk.table = TRUE)

从图5看到，曲线是反过来的，但是仔细看纵坐标不是0-1，原来”cumhaz” 绘制的是随时间推移事件发生的累计风险(HR)。

其实，仔细参阅ggsurvplot的参数说明，这个问题就迎刃而解了：

- "cumhaz" plots the cumulative hazard function (f(y) = -log(y)), 绘制累计风险图
- "event" plots cumulative events (f(y) = 1-y), 绘制1-S(t)曲线
- "pct" for survival probability in percentage，绘制百分比生存概率图

定义fun = "event"，输出置信区间，接下来调整一下图的格式、配色，就可以得到一张漂亮的事件发生的累计概率曲线图啦。

ggsurvplot(survfit(Surv(Os,Status)~Group,data=data1),
           fun = "event", # l绘制1-S(t)曲线
           conf.int = TRUE,  # 输出曲线置信区间
           conf.int.style="ribbon", # 可定义的类别包括("ribbon", "step")
           conf.int.alpha=0.45, # 设置置信区间透明区
           surv.plot.height= 0.7, # 生存图高度
           main = "Survival curve",
           legend.title = "Group",
           legend.labs = c("Treat", "Control"),
           legend=c(0.8,0.7),
           xlab="Days",
           ylab="Cumulative Event Rate (%)", # 或者直接写“1-S(t)”
           xlim = c(0,60), ylim = c(0,1),
           risk.table = TRUE,
           tables.height = 0.15, # 生存曲线图下所有生存表的高度，数值0-1之间
           tables.theme = theme_cleantable(),
           font.main = c(14, "bold", "darkblue"),
           font.x = c(14, "plain", "black"),
           break.x.by = 10, break.y.by = 0.1,
           font.y = c(14, "plain", "black"),
           font.tickslab = c(12, "plain", "black"),
           risk.table.col="strata",
           risk.table.height=0.2,
           palette = "npg", # 可选调色板有 "grey","npg","aaas","lancet","jco","ucscgb","uchicago","simpsons"和"rickandmorty".
           ggtheme = theme_bw() # Change ggplot2 theme
          )