受限平均存活时间（RMST）及其 R 语言实现。
前些天我的学徒写了教程：人人都可以学会生存分析（学徒数据挖掘） 吸引到了读者：武汉大学金文意，他希望可以分享一下生存分析的新玩法：
受限平均生存时间
大家好，我是武汉大学金文意。很荣幸受曾老师邀请，与大家简短地分享一下我对受限平均生存时间的一点薄见。水平有限，有不足之处敬请指正！
背景
在诸如JAMA oncology等*期刊中，我们经常会看到如图1所示的Restricted mean survival time（RMS time），即受限平均生存时间1。基于RMS time的组间评价方法既往多见于前瞻性的临床对照实验中，但是近年来，RMS time在预后模型（prognostic signature）领域亦逐渐开始普及。
单从与生物信息学息息相关的prognostic signature领域来说，假如你构建了一个prognostic signature，该signature可将患者分成不同风险组（例如高、中、低风险）。传统的做法是直接绘制K-M生存曲线并log-rank计算组间差异，如果p-value小于0.05即可判断组间的生存率有显著差异。
但是，仅仅这样就足够了吗？log-rank test比较的是存在删失的情况下组间的生存时间分布2。对于更精确的医学研究，我们更想知道的是，截止至某一年的时间点，组间的生存差异。这个时候，就要引出我们的RMS time了。
图1
定义
RMS time, 其描述的是纳入研究的患者在随访t时间段内的平均生存时间，可以估计为生存曲线下直至预定时间点的面积，因此可以考虑该时间范围之前的所有生存信息3。这恰恰是其相对于传统生存分析的优势所在。
R语言实现
survRM2 package常用来计算RMS time，方便起见，使用包中的示例数据进行展示。
require(survRM2)
data =rmst2.sample.data() #input data
我们看一下数据结构，我们只需要三列数据（生存时间、生存状态、及你的依据处理因素/风险因素的分组arm）
head(data[,1:3])## time status arm## 1 1.095140 1 1## 2 12.320329 0 1## 3 2.770705 1 1## 4 5.270363 1 1## 5 4.117728 0 0## 6 6.852841 1 0简单调用函数进行计算,依序输入，其中，tau指的是你所设置的时间点。函数的其他参数的具体说明可自行查看R中该包的说明手册。
RMS <- rmst2(data$time, data$status, data$arm, tau=10)print(RMS)## The truncation time: tau = 10 was specified.## ## Restricted Mean Survival Time (RMST) by arm## Est. se lower .95 upper .95## RMST (arm=1) 7.146 0.283 6.592 7.701## RMST (arm=0) 7.283 0.295 6.704 7.863###### Restricted Mean Time Lost (RMTL) by arm## Est. se lower .95 upper .95## RMTL (arm=1) 2.854 0.283 2.299 3.408## RMTL (arm=0) 2.717 0.295 2.137 3.296###### Between-group contrast## Est. lower .95 upper .95 p## RMST (arm=1)-(arm=0) -0.137 -0.939 0.665 0.738## RMST (arm=1)/(arm=0) 0.981 0.878 1.096 0.738## RMTL (arm=1)/(arm=0) 1.050 0.787 1.402 0.738可以看到给出的信息非常详细并给出了组间的显著性指标，一般来说，文章中只需要给出RMS time即可。
可视化的话，可以直接使用plot函数，也可依据你的个人需求编制函数进行定制,将其组合成一幅图等。我们简单的plot一下：
plot(RMS, xlab="Years", ylab="Probability")其中，红色部分表示RMS time，黄色部分表示RMSL（Restricted Mean Time Lost）。
总体来说，RMS time的分析过程较为简单。然而，RMS time这项指标，如同时间依赖的ROC曲线一样，恰恰是佐证你研究的一剂强心剂。
参考文献
1. Li B, Cui Y, Diehn M, Li R. Development and Validation of an Individualized Immune Prognostic Signature in Early-Stage Nonsquamous Non-Small Cell Lung Cancer. JAMA Oncol. 2017;3(11):1529-1537. doi:10.1001/jamaoncol.2017.1609
2. http://www.lifescienceglobal.com/pms/index.php/ijsmr/article/view/674
3. Wang Z, Wu H, Xie L, et al. Correlation of Milestone Restricted Mean Survival Time Ratio With Overall Survival Hazard Ratio in Randomized Clinical Trials of Immune Checkpoint Inhibitors: A Systematic Review and Meta-analysis. JAMA Netw Open. 2019;2(5):e193433. doi:10.1001/jamanetworkopen.2019.3433
写在文末
如果你确实想把生存分析学会，学好，我在生信技能树多次分享过生存分析的细节也值得你认真读完；
- 基因表达量高低分组的cox和连续变量cox回归计算的HR值差异太大?
- 学徒作业-两个基因突变联合看生存效应
- TCGA数据库里面你的基因生存分析不显著那就TMA吧
- 批量cox生存分析结果也可以火山图可视化
- 既然可以看感兴趣基因的生存情况，当然就可以批量做完全部基因的生存分析
- 多测试几个数据集生存效应应该是可以找到统计学显著的！
- 我不相信kmplot这个网页工具的结果(生存分析免费做)
- 为什么不用TCGA数据库来看感兴趣基因的生存情况
- 200块的代码我的学徒免费送给你，GSVA和生存分析
- 集思广益-生存分析可以随心所欲根据表达量分组吗
- 生存分析时间点问题
- 寻找生存分析的最佳基因表达分组阈值
- apply家族函数和for循环还是有区别的（批量生存分析出图bug）
- TCGA数据库生存分析的网页工具哪家强
- KM生存曲线经logRNA检验后也可以计算HR值