R言語-生存分析

生存分析に関する知識は、自分で「アルゴリズムとモデル」類のブログで勉強してください.(http://blog.csdn.net/xiaohukun/article/details/77679134)
一、パッケージをダウンロードしてロードする
生存分析にはsurvivalパッケージを使用しています

install.packages("survival")   #  survival 
library(survival)      #  survival 

二、データの準備
survivalパケットに付属する「pbc」データセットを例(418*20)に記録されたtimeは元の時間にすぎず、statusに基づいて生存時間のタイプを判断し、対応するフォーマットに変換する必要がある.この仕事はSurv(time,event)が完成し、生存対象に戻る.

> Sur_Obj 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にベクトルが1つしかない場合は、x~1の形式に書きます.

#      
> model1) 
> 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=" ")を直接通過すると生存曲線を描くことができます.
四、単因分析
ここでは主に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 

元の仮定は2種類のデータの平均値に明らかな差がなく、P値は元の仮定を受信する確率であり、明らかにここで>>0.05であるため、元の仮定を拒否することはできない.結論:trtの異なる患者の生存時間には明らかな差はなかった.
五、多要素分析
逐次回帰方式を用いて,影響因子がモデルに入るか否かを一歩一歩決定し,コヒーレント変数のPH仮定検証を行った.coxph(formula)関数とcoxをそれぞれ用いる.zph(fit)関数.

#     ，            ，P         ，coef   
> coxmodel 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 zphmodel
             rho  chisq     p
pbc$bili 0.00717 0.0064 0.936

描画することもできますβ 図の方法は係数の変化傾向を見て、変化が大きくなければPH仮定にも合致すると考えられる.plot(zphmodel)
このようにして,生存分析に関わる基本分析が完了した.