生存分析 R语言（一）——KM estimator

题目

Data contains drug users who have been given residential treatment for drug abuse. Consider only four columns：age, time, treat, censor. Treat and censor are indicators.
The main task is all about plot the KM estimator and test between several groups. Here focus on the code only.

packages的安装与加载

安装R包，只需在第一次时安装即可，后续使用无须再次安装

install.packages("survminer")
install.packages("survival")

加载R包，每次使用均需加载

library（survminer)
library(survival)

数据的导入

方案一

将数据直接放入工作空间中

getwd()#可直接得到当前工作空间地址
setwd()#可自行设置工作空间地址
read.csv("file.csv",header=TRUE)#header默认为真，即第一行为变量名

方案二

mydata<-read.csv(file.choose())#回车手动选择文件

Kaplan-Meier estimator

拟合曲线

attach(mydata)#绑定数据集，当环境中多个数据时很必要
Surv(time,censor)#创建生存对象

加号为censoring data

  [1]   23   167   186   232   550+   71   398     3   559    61    75   548+  256 
 [14]   51    83   597+  190   208     8   130   144    20    74   119   175   567+
 [27]  546+   40    76   461    12   224   211   597+  619+  124    11   211   353 
 [40]  551+   89    65    37   104   190    55   490    56    35   603+  362   116 
 [53]    7   119   217   575+  134    54   233   540+  238   339    90   242   116 
 [66]   68   228   655+  122    26   569+  137   203   364    59   198   166    23 
 [79]   85    23   271    99   180    99   549+  181   103   203   519   541+   11 
 [92]   11   119    13   307    73   373   164   436   593+  292    54    45    53 
[105]   51   122   121   556+  554+  245   159   268   123    74    86   538+  131 
[118]  538+  434   106   546+  533+   12   537+  396   248    22    84    37    79 
[131]   45   121   168   567+   93   630+  194   317   181     9    14   539+   14 
[144]  403    28   510+   81   547+  226     9    58   140    72   546+  548+  129 
[157]  231    86   196    19    21   624+   37    15   127   105    35   267    59 
[170]  115   222   282   120   148   220   180   102   220   174   548+  231   348 
[183]  386   198   123   548+  182   567+   35    87    53   227    47     5   763+
[196]  374   148    94    65   634+   89    90   551+  437   136   290   159   354 
[209]   96   659   260   105    81   170    78   138    71   139   280    32   289 
[222]   57   559+  122   156   166   805+  107   169    29   171   248   118   164 
[235]   51   104   537+    6    92   501+   87   499    78    44   361   102    36 
[248]   94    76   133   184    31   210    54    46   123   720+  250   106   226 
[261]   26   110   279   541+  533+  545+  565+  147    44    69   548+  245    45 
[274]  177   130   551+  110    56    14   354    51   232    68   144   598+  449 
[287]   79   125    90    87   633+  127   224   188    56   285   156   578+  313 
[300]   25    36     3   155   122   471   539+  186   285    52   227    97    10 
[313]  431   224   360   551+   26   207   659+   87     5   151   220   167   762+
[326]    2    70    84   162   567+  477   299    21   562+  325   540+  115   545+
[339]  243   199   170   100   130    85   540+  157   287    30   399    80   102 
[352]  274   106   192   510+   85    78   115   562+   33   384   587+  541+  228 
[365]  560+  244    74    50    79    91   153    15   156   259   560+  494   276 
[378]   32   142   567+  560+  144    69   136   543+   24    73    10   306    29 
[391]  557+  502   540+   44    79   549+  547+    5   225    87   246     6   321 
[404] 1172+  734+  285    87   544+  253   442   625+   99   536+  621+   79   569+
[417]  575+  176   576+  150   408    81   191     7   540+  654+  120   379   184 
[430]  441    44   226   297    38   343   491    67    70    46   568+  194   239 
[443]  328    92   164     4    15    98   546+   58    95    65   459   193   102 
[456]  349   192   503+  547+  120   567+  226   458   554+   51   259    43    31 
[469]   30     4     4    90   547+  354   115   142   307    86    23    26     5 
[482]  121   183    57    19    38    84    21   568+  124   216     7   115    32 
[495]  161   200    47    48    60    18   212   231    96    83   586+  380    92 
[508]   57   152   209    14   187    81   332   101   547+  168   243   154   255 
[521]  401   142    10   188   459   523+   95   273   218   132     4   115   113 
[534]   91   143    20   369   107    59   255   159   168   210   168    94   106 
[547]  552+  242   113   555+   23    62   106   176    98   555+  290   503+  548+
[560]  233   119   542+   34   170   215    46   296   260   655+   67   274   537+
[573]  516   268    69   355   184   545+    7   544+  560+   36    47   589+  274 
[586]   83   658+  175   555+  175   568   538+  306     6   582+  220   103   708+
[599]  100   324   300    22    80   237   161   159   541+  162    65   143   556+
[612]  204   538+  310   552+  377   111   564+   24   189   547+  126   553+    8 
[625]   17   540+    3    10 

对象进行拟合

fit<-survfit(Surv(time,censor)~treat,data=mydata)#拟合生存曲线
fit#可查看拟合曲线信息

Call: survfit(formula = Surv(time, censor) ~ treat, data = mydata)

          n events median 0.95LCL 0.95UCL
treat=0 320    265    132     115     156
treat=1 308    243    192     175     226

可以看到n是每组内的人数，events为事件（drug abuse)发生次数，后面分别是中位数生存时间、95%置信区间的下界限和上界限。

绘制拟合曲线

ggsurvplot(fit,data=mydata,
           surv.median.line="hv",
           conf.int=TRUE)
           #分别增加了中位数和置信区间
           #palette可自定义调色板，默认为hue.

假设检验比较两组生存曲线的差异

survdiff(Surv(time,censor)~treat,data=mydata,rho=0)

rho代表了假设检验中统计量的weight。
rho=0 Log-Rank test
rho=1 Generalized Wilcoxon test
rho=1/2 Tarone-Ware test

Call:
survdiff(formula = Surv(time, censor) ~ treat, data = mydata)

          N Observed Expected (O-E)^2/E (O-E)^2/V
treat=0 320      265      236      3.61       6.8
treat=1 308      243      272      3.13       6.8

 Chisq= 6.8  on 1 degrees of freedom, p= 0.009