R语言-生存分析与结果的图像处理

数据准备：

library("survival")
library("survminer")
data("lung")

调用“lung”数据集，使用head()命令调查前6行，得到以下结果：

 inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
1    3  306      2  74   1       1       90       100     1175      NA
2    3  455      2  68   1       0       90        90     1225      15
3    3 1010      1  56   1       0       90        90       NA      15
4    5  210      2  57   1       1       90        60     1150      11
5    1  883      2  60   1       0      100        90       NA       0
6   12 1022      1  74   1       1       50        80      513       0

time是指生存时间，status是二值型的生存状态，其中“1”代表存活，“2”代表已经死亡。

KM非参数法：

调用survival包中的survfit来进行生存分析，得到每个时间段中的生存率。

fit <- survfit(Surv(time, status) ~ 1, data = lung)

直接查看fit的信息可以看到

> fit
Call: survfit(formula = Surv(time, status) ~ 1, data = lung)n  events  median 0.95LCL 0.95UCL 228     165     310     285     363

n:代表总的样本数为228个
events:代表status为2的样本数（死亡样本数）
0.95LCL:95% 置信区间的下界
0.95UCL:95%置信区间的上界
median:代表中位数

Log-Rank检验：

针对使用KM方法的生存分析，可以使用survdiff()来比较两组或者多组之间的的生存曲线，调查它们生存率和标准误之间的差异，输出p值。

> survdiff(Surv(time, status) ~ ph.karno, data = lung)
Call:
survdiff(formula = Surv(time, status) ~ ph.karno, data = lung)n=227, 1 observation deleted due to missingness.N Observed Expected (O-E)^2/E (O-E)^2/V
ph.karno=50   6        5     5.71    0.0874    0.0935
ph.karno=60  19       16     9.64    4.1905    4.5095
ph.karno=70  32       29    20.53    3.4943    4.0275
ph.karno=80  67       47    43.30    0.3161    0.4390
ph.karno=90  74       49    58.85    1.6489    2.5910
ph.karno=100 29       18    25.97    2.4454    2.9262Chisq= 12.3  on 5 degrees of freedom, p= 0.03

COX半参数法：

使用coxph()函数来进行生存分析：

cox_single <- coxph(Surv(time, status) ~ sex, data = lung)    #单变量
cox_muti <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data =  lung)#多变量

使用summary函数可以查看返回当结果和p值：

> summary(cox_single)
Call:
coxph(formula = Surv(time, status) ~ 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 ‘ ’ 1exp(coef) exp(-coef) lower .95 upper .95
sex     0.588      1.701    0.4237     0.816Concordance= 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

coxph()会给出Likelihood ratio test，Wald test，Score (logrank) test三种检验的结果，当样本数量足够大时，三者的值相似，一般情况下Likelihood ratio test的值更加准确。

相似比检验：

针对cox模型，R里面的cox.zph函数可以对其进行检验：

R<-cox.zph(cox.muti)> Rchisq df    p
age     0.188  1 0.66
sex     2.305  1 0.13
ph.ecog 2.054  1 0.15
GLOBAL  4.464  3 0.22

其原假设为符合PH假定，代表假定HR值不随时间发生变化，满足比例风险假定。在结果中所有变量的p值均>0.1，无法拒绝原假设。

生存曲线和诊断图：

1，plot（）
对于KM方法而言可以直接使用plot()命令来绘制生存曲线:

plot(fit,main="survival-KM",xlab="time",ylab="rate")

2，ggsurvplot
使用ggsurvplot（）函数对KM的结果绘制生存曲线：

ggsurvplot(fit,pval = FALSE, conf.int = TRUE, #是否显示p值/显示风险区域risk.table = TRUE, # 将风险表显示在生存曲线下面risk.table.col = "strata", linetype = "strata", surv.median.line = "hv", # 中位数的标记可选（hv,v,h）ggtheme = theme_bw(), # 改变ggplot2的样式palette = "green"    #设置曲线颜色
)

3，多组生存曲线
在绘制多组间生存曲线时，pval=TRUE可以输出log-rank的p值，使用pval.method=T可以显示计算p值的方法；

调查ph.ecog组间生存曲线的差异：

fit <- survfit(Surv(time, status) ~ ph.ecog, data = lung)
ggsurvplot(fit,data=lung,pval=T,risk.table = T,risk.table.col="red",pval.method=T)

%其他参数：

risk.table = TRUE    #风险表是否显示
risk.table.col="red" #风险表颜色
surv.median.line = "hv"#中位数标记的方法（hv,h,v）

4，COX模型的Hazard ratio图：
使用ggforest()函数进行绘制，可以得到log-rank的p值，AIC值
5，Schoenfeld individual检验诊断图：
又称“Schoenfeld Residuals Test”，用于调查残差与时间的独立性，从而检验Cox模型中的比例风险假设。相当于是cox.zph()结果的可视化表现。P值不显著，说明符合原假设

R<-cox.zph(cox_muti)  #cox多变量检验
ggcoxzph(R,resid=T,se=T,point.col = "red",point.size = 1)

参考文献：

https://www.mbaskool.com/business-concepts/statistics/8766-schoenfeld-residuals-test.html
https://blog.csdn.net/weixin_43700050/article/details/100938707
https://blog.csdn.net/xiaohukun/article/details/78019726
https://blog.csdn.net/jaen_tail/article/details/79081954

*刚刚接触生存分析这一块，还有很多东西需要去学习，望大家一起讨论一起学习。Fight!