R语言统计与绘图：生存曲线的多重比较-经管之家官网！

扫码加入统计交流群

来源：R语言统计与绘图 公众号，做学习分享用。

前面我们学习了基础绘图包和ggplot2包绘制生存曲线，还有生存率的比较，今天来学习下生存曲线的两两比较。

survdiff()函数可以检验两条或者多条生存曲线之间是否有差异，但这个比较的是总体间是否存在差异，不能比较生存曲线之间两两是否存在差异。

在R中使用survminer包的pairwise_survdiff()函数用来进行生存曲线之间的多重比较。

1. 加载数据集

library(survival) # 加载含内置数据集的包
data(colon) # 加载数据集
View(colon) # 预览数据集
str(colon) # 查看数据集结构

数据集介绍见这篇文章——常用内置数据集介绍。

从数据集介绍中我们知道rx为三分类变量，extent为四分类变量，status为删失状态，time为生存时间。

2. 生存曲线的多重比较(一)

根据rx变量分组绘制生存曲线。

library(survminer) # 加载包
restest <- pairwise_survdiff(Surv(time, status) ~ rx,
                             data = colon)
restest

输出：
      Pairwise comparisons using Log-Rank test 
data:  colon and rx 
        Obs     Lev    
Lev     0.78    -      
Lev+5FU 3.3e-07 9.5e-07
P value adjustment method: BH 

结果解释：
rx变量共分为三组：观察组(Obs)，治疗组1(Lev)和治疗组2(Lev+5FU)。
从结果可知，Obs组 vs Lev组 p-value=0.78，无统计学意义，两组生存率没有差异。同理可知，Obs组 vs Lev+5FU组和Lev组 vs Lev+5FU组的生存率都有统计学差异。

3. 生存曲线的多重比较(二)

根据extent变量分组绘制生存曲线。

在colon数据集，extent为四分类变量，我们先将其转换为因子，并添加标签。

colon$extent <- factor(colon$extent, 
                       levels = c(1,2,3,4),
                       label = c("submucosa","muscle","serosa","contiguous structures"))
is.factor(colon$extent) # 查看变量是否为因子
[1] TRUE
levels(colon$extent) # 查看因子水平
[1] "submucosa"  "muscle"   "serosa"               
[4] "contiguous structures"

# library(survminer) # 加载包
restest <- pairwise_survdiff(Surv(time, status) ~ extent,
                             data = colon)
restest

输出：
      Pairwise comparisons using Log-Rank test 
data:  colon and extent 
                      submucosa muscle  serosa 
muscle                0.16535   -       -      
serosa                0.00064   8.1e-07 -      
contiguous structures 2.8e-06   7.1e-11 0.00013
P value adjustment method: BH

结果解释：同rx。

4. P值用符号"*"表示

此函数可以使用符号“*”来表示p值范围。

symnum(restest$p.value, cutpoints = c(0, 0.0001, 0.001, 0.01, 0.05, 0.1, 1),
       symbols = c("****", "***", "**", "*", "+", " "),
       abbr.colnames = FALSE, na = "")

输出：
                      submucosa muscle serosa
muscle                                       
serosa                ***       ****         
contiguous structures ****      ****   ***   
attr(,"legend")
[1] 0 ‘****’ 1e-04 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1 \t    ## NA: ‘’

说明：
p-value在0-0.0001之间用“四颗*”表示；
在0.0001-0.001之间用“三颗*”表示；
在0.001-0.01之间用“两颗*”表示；
在0.01-0.05之间用“一颗*”表示；
在0.05-0.1之间用“+”表示；
在0.1-1之间用“空值NA”表示

5. pairwise_survdiff()函数

pairwise_survdiff(formula, data, p.adjust.method = "BH", na.action,
                  rho = 0)

formula：与其他生存模型一样的表达式，表示为Surv(time, status) ~ predictors。  
data：数据框，公式中变量的来源。 

p.adjust.method：调整p值的方法（参阅p.adjust）。
允许的值包括"holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none".
如果不想调整p值(不建议)，请使用p.adjust.method = "none"。  

na.action：缺失值处理函数。  
rho：控制检验类型的参数，允许的值有0（用于log-rank检验）和1（用于peto & peto检验）。  

参考资料：
pairwise_survdiff()函数帮助文件

End

「经管之家」APP：经管人学习、答疑、交友，就上经管之家！
免流量费下载资料----在经管之家app可以下载论坛上的所有资源，并且不额外收取下载高峰期的论坛币。
涵盖所有经管领域的优秀内容----覆盖经济、管理、金融投资、计量统计、数据分析、国贸、财会等专业的学习宝库，各类资料应有尽有。
来自五湖四海的经管达人----已经有上千万的经管人来到这里，你可以找到任何学科方向、有共同话题的朋友。
经管之家（原人大经济论坛），跨越高校的围墙，带你走进经管知识的新世界。
扫描下方二维码下载并注册APP