R语言-生存分析

淡淡的烟草味﹌ 2022-06-08 06:46 345阅读 0赞

有关生存分析的相关知识，请自行前往“算法与模型”类博客进行学习。  
([http://blog.csdn.net/xiaohukun/article/details/77679134][http_blog.csdn.net_xiaohukun_article_details_77679134])

#### 一、下载并加载包 ####

进行生存分析使用的是**survival** 包

install.packages("survival")   #安装survival包
    library(survival)      #加载survival包

#### 二、准备数据 ####

使用survival包自带的“pbc”数据集为例（418\*20）  
![这里写图片描述][SouthEast]  
这里记录的time只是原始时间，还需要根据status判断生存时间的类型，转化成相应的格式。这个工作由**Surv(time,event)**完成,并返回生存对象

> Sur_Obj<-Surv(pbc$time,pbc$status)
    > Sur_Obj
    [1]  400  4500? 1012  1925  1504+ 2503  1832? 2466  2400 
     [10]   51  3762   304  3577? 1217  3584  3672?  769   131 ……

#### 三、绘制生存曲线 ####

函数**surfit(formula)**，根据填入的formula按照默认K-M算法进行生存曲线拟合。如果填入的formula只有一个向量，则写成x~1的形式。

#获得曲线模型
    > model<-survfit(Sur_Obj~1) 
    > summary(model)
    Call: survfit(formula = Sur_Obj ~ 1)
    
    232 observations deleted due to missingness 
     time n.risk n.event survival std.err lower 95% CI upper 95% CI
       41    186       2  0.98925 0.00756      0.97454       1.0000
       43    184       1  0.98387 0.00924      0.96593       1.0000
       51    183       1  0.97849 0.01064      0.95787       0.9996
       71    182       1  0.97312 0.01186      0.95015       0.9966
       77    181       1  0.96774 0.01296      0.94268       0.9935
       94    180       1  0.96237 0.01395      0.93540       0.9901
      110    179       1  0.95699 0.01488      0.92827       0.9866
      111    178       1  0.95161 0.01573      0.92127       0.9830
      130    177       1  0.94624 0.01654      0.91437       0.9792
      ……
    #其中lower 95% CI 和 upper 95% CI为95%置信度时生存率的上下区间
    
    > model
    Call: survfit(formula = Sur_Obj ~ 1)
    
       232 observations deleted due to missingness 
          n  events  median 0.95LCL 0.95UCL 
        186     161    1217    1077    1492 
    #其中median 即为中位生存期,0.95LCL,0.95UCL为预测上下限

直接通过`plot(model,ylab = "生存率",xlab="天")`便可以绘制生存曲线。  
![这里写图片描述][SouthEast 1]

#### 四、单因素分析 ####

这里主要是指用**survdiff(formula)**函数进行log-rank检验。

> survdiff(Sur_Obj~pbc$trt)   #trt是分组条件
    Call:
    survdiff(formula = Sur_Obj ~ pbc$trt)
    
    n=144, 274 observations deleted due to missingness.
    
               N Observed Expected (O-E)^2/E (O-E)^2/V
    pbc$trt=1 75       65     65.5   0.00441   0.00956
    pbc$trt=2 69       60     59.5   0.00486   0.00956
    
     Chisq= 0  on 1 degrees of freedom, p= 0.922

原假设为两类数据均值无明显差异，P值为接收原假设的概率，显然此处>>0.05，因此不能拒绝原假设。结论：trt不同的患者生存时间并无明显差异。

#### 五、多因素分析 ####

采用逐步回归的方式一步一步确定影响因素是否进入模型，并进行协变量的PH假定验证。分别使用到**coxph(formula)**函数和**cox.zph(fit)**函数。

#相似比检验，原假设为加入新因素无差别，P值看是否为影响因素，coef为系数
    > coxmodel<-coxph(Sur_Obj~pbc$bili)
    > coxmodel
    Call:
    coxph(formula = Sur_Obj ~ pbc$bili)
    
              coef exp(coef) se(coef)    z       p
    pbc$bili 0.078     1.081    0.013 6.01 1.8e-09
    
    Likelihood ratio test=28.6  on 1 df, p=9.09e-08
    n= 186, number of events= 161 
       (232 observations deleted due to missingness)
    
    #原假设为符合PH假定 
    > zphmodel<-cox.zph(coxmodel)
    > zphmodel
                 rho  chisq     p
    pbc$bili 0.00717 0.0064 0.936

也可以通过绘制β图的方法看系数的变化趋势，如变化不大则也认为符合PH假定。  
`plot(zphmodel)`  
![这里写图片描述][SouthEast 2]

如此，生存分析中涉及的基本分析就完成了。

[http_blog.csdn.net_xiaohukun_article_details_77679134]: http://blog.csdn.net/xiaohukun/article/details/77679134
[SouthEast]: /images/20220608/83da42f9491244c88bdccd7cf78b5c7d.png
[SouthEast 1]: /images/20220608/bb475152efad48efbd9d63bc8560deae.png
[SouthEast 2]: /images/20220608/2a0fc216619d43dea68d6636b967c9be.png