生存分析是临床常用统计方法，一旦和时间扯上关系，分析就变得复杂多了，此时不再是单一的因变量，还需要考虑时间给因变量和自变量带来的各种影响。
本次主要演示R语言做生存分析的一些方法。比如寿命表、K-M曲线、logrank检验。后续还会给大家介绍Cox回归、时依系数和时依协变量的Cox回归、生存曲线的可视化 等内容。
本推文不涉及理论，只有实操，想要了解生存分析的理论的请自行学习。不涉及理论，并不代表理论不重要，在以后的机器学习和临床预测模型的相关推文中，会经常用到这些理论，建议大家学习一下。
生存过程的描述 library (survival)library (survminer)使用survival包中的lung数据集用于演示，这是一份关于肺癌患者的生存数据。time是生存时间，以天为单位，status是生存状态，1代表删失，2代表死亡。但是一般在生存分析中我们喜欢用1代表死亡，用0代表删失，所以我们更改一下（其实不改也可以，你记住就行）。
df <- lung2 ,1 ,0 )## 'data.frame': 228 obs. of 10 variables: ## $ inst : num 3 3 3 5 1 12 7 11 1 7 ... ## $ time : num 306 455 1010 210 883 ... ## $ status : num 1 1 0 1 1 0 1 1 1 1 ... ## $ age : num 74 68 56 57 60 74 68 71 53 61 ... ## $ sex : num 1 1 1 1 1 1 2 2 1 1 ... ## $ ph.ecog : num 1 0 0 1 0 1 2 2 1 2 ... ## $ ph.karno : num 90 90 90 90 100 50 70 60 70 70 ... ## $ pat.karno: num 100 90 90 60 90 80 60 80 80 70 ... ## $ meal.cal : num 1175 1225 NA 1150 NA ... ## $ wt.loss : num NA 15 15 11 0 0 10 1 16 34 ... 首先把生存时间和生存状态用Surv()放到一起，可以看到有+的就是截尾数据。
Surv(time = lung$time, event = lung$status)## [1] 306 455 1010+ 210 883 1022+ 310 361 218 166 170 654 ## [13] 728 71 567 144 613 707 61 88 301 81 624 371 ## [25] 394 520 574 118 390 12 473 26 533 107 53 122 ## [37] 814 965+ 93 731 460 153 433 145 583 95 303 519 ## [49] 643 765 735 189 53 246 689 65 5 132 687 345 ## [61] 444 223 175 60 163 65 208 821+ 428 230 840+ 305 ## [73] 11 132 226 426 705 363 11 176 791 95 196+ 167 ## [85] 806+ 284 641 147 740+ 163 655 239 88 245 588+ 30 ## [97] 179 310 477 166 559+ 450 364 107 177 156 529+ 11 ## [109] 429 351 15 181 283 201 524 13 212 524 288 363 ## [121] 442 199 550 54 558 207 92 60 551+ 543+ 293 202 ## [133] 353 511+ 267 511+ 371 387 457 337 201 404+ 222 62 ## [145] 458+ 356+ 353 163 31 340 229 444+ 315+ 182 156 329 ## [157] 364+ 291 179 376+ 384+ 268 292+ 142 413+ 266+ 194 320 ## [169] 181 285 301+ 348 197 382+ 303+ 296+ 180 186 145 269+ ## [181] 300+ 284+ 350 272+ 292+ 332+ 285 259+ 110 286 270 81 ## [193] 131 225+ 269 225+ 243+ 279+ 276+ 135 79 59 240+ 202+ ## [205] 235+ 105 224+ 239 237+ 173+ 252+ 221+ 185+ 92+ 13 222+ ## [217] 192+ 183 211+ 175+ 197+ 203+ 116 188+ 191+ 105+ 174+ 177+ 如果只是想要描述一下这份生存数据，可以使用寿命表法或者K-M曲线，在R中可以通过survfit()实现。
# 构建生存曲线 1 , data = df)# 寿命表，surv_summary比默认的summary()更好 ## time n.risk n.event n.censor surv std.err upper lower ## 1 5 228 1 0 0.99561404 0.004395615 1.0000000 0.98707342 ## 2 11 227 3 0 0.98245614 0.008849904 0.9996460 0.96556190 ## 3 12 224 1 0 0.97807018 0.009916654 0.9972662 0.95924368 ## 4 13 223 2 0 0.96929825 0.011786516 0.9919508 0.94716300 ## 5 15 221 1 0 0.96491228 0.012628921 0.9890941 0.94132171 ## 6 26 220 1 0 0.96052632 0.013425540 0.9861367 0.93558107 ## 7 30 219 1 0 0.95614035 0.014184183 0.9830945 0.92992527 ## 8 31 218 1 0 0.95175439 0.014910735 0.9799794 0.92434234 ## 9 53 217 2 0 0.94298246 0.016284897 0.9735659 0.91335978 ## 10 54 215 1 0 0.93859649 0.016939076 0.9702809 0.90794671 ## 11 59 214 1 0 0.93421053 0.017574720 0.9669508 0.90257880 ## 12 60 213 2 0 0.92543860 0.018798180 0.9601711 0.89196244 ## 13 61 211 1 0 0.92105263 0.019389168 0.9567281 0.88670745 ## 14 62 210 1 0 0.91666667 0.019968077 0.9532533 0.88148430 ## 15 65 209 2 0 0.90789474 0.021093908 0.9462168 0.87112471 ## 16 71 207 1 0 0.90350877 0.021642644 0.9426590 0.86598451 ## 17 79 206 1 0 0.89912281 0.022182963 0.9390770 0.86086855 ## 18 81 205 2 0 0.89035088 0.023240987 0.9318456 0.85070391 ## 19 88 203 2 0 0.88157895 0.024272607 0.9245323 0.84062118 ## 20 92 201 1 1 0.87719298 0.024779731 0.9208475 0.83560803 ## 省略一部分结果。。。 ## 178 791 9 1 0 0.07831533 0.314346559 0.1450170 0.04229358 ## 179 806 8 0 1 0.07831533 0.314346559 0.1450170 0.04229358 ## 180 814 7 1 0 0.06712742 0.350176075 0.1333431 0.03379322 ## 181 821 6 0 1 0.06712742 0.350176075 0.1333431 0.03379322 ## 182 840 5 0 1 0.06712742 0.350176075 0.1333431 0.03379322 ## 183 883 4 1 0 0.05034557 0.453824434 0.1225342 0.02068546 ## 184 965 3 0 1 0.05034557 0.453824434 0.1225342 0.02068546 ## 185 1010 2 0 1 0.05034557 0.453824434 0.1225342 0.02068546 ## 186 1022 1 0 1 0.05034557 0.453824434 0.1225342 0.02068546 画出生存曲线，横坐标是生存时间，纵坐标是生存率。
ggsurvplot(fit,TRUE , # 可信区间 'blue' , # 更改配色 "hv" , # 中位生存时间 # 更改主题 生存过程的比较 如果通过某个变量把数据分为多组，然后检验不同组别之间的生存时间（生存曲线）有无差别，则可以通过logrank检验或者breslow检验。
在R语言中通过survdiff()实现logrank检验。
fit <- survdiff(Surv(time, status) ~ sex, data = df)## Call: ## survdiff(formula = Surv(time, status) ~ sex, data = df) ## ## N Observed Expected (O-E)^2/E (O-E)^2/V ## sex=1 138 112 91.6 4.55 10.3 ## sex=2 90 53 73.4 5.68 10.3 ## ## Chisq= 10.3 on 1 degrees of freedom, p= 0.001 可以用神包broom提取数据：
broom::tidy(fit)## # A tibble: 2 × 4 ## sex N obs exp ## <chr> <dbl> <dbl> <dbl> ## 1 1 138 112 91.6 ## 2 2 90 53 73.4 broom::glance(fit)## # A tibble: 1 × 3 ## statistic df p.value ## <dbl> <dbl> <dbl> ## 1 10.3 1 0.00131 对于不同组别之间生存曲线的检验，也可以通过K-M图示的方法：
fit.logrank <- survfit(Surv(time, status) ~ sex, data = df)# 可以查看寿命表 ## time n.risk n.event n.censor surv std.err upper lower ## 1 11 138 3 0 0.97826087 0.01268978 1.0000000 0.95423012 ## 2 12 135 1 0 0.97101449 0.01470747 0.9994124 0.94342350 ## 3 13 134 2 0 0.95652174 0.01814885 0.9911586 0.92309525 ## 4 15 132 1 0 0.94927536 0.01967768 0.9866017 0.91336116 ## 5 26 131 1 0 0.94202899 0.02111708 0.9818365 0.90383547 ## 6 30 130 1 0 0.93478261 0.02248469 0.9768989 0.89448204 ## 7 31 129 1 0 0.92753623 0.02379334 0.9718155 0.88527450 ## 8 53 128 2 0 0.91304348 0.02627035 0.9612864 0.86722163 ## 9 54 126 1 0 0.90579710 0.02745220 0.9558688 0.85834834 ## 10 59 125 1 0 0.89855072 0.02860313 0.9503632 0.84956295 ## 11 60 124 1 0 0.89130435 0.02972717 0.9447781 0.84085714 ## 12 65 123 2 0 0.87681159 0.03190746 0.9333961 0.82365740 ## 13 71 121 1 0 0.86956522 0.03296902 0.9276100 0.81515252 ## 14 81 120 1 0 0.86231884 0.03401448 0.9217668 0.80670491 ## 15 88 119 2 0 0.84782609 0.03606427 0.9099232 0.78996675 ## 16 92 117 1 0 0.84057971 0.03707173 0.9039292 0.78166991 ## 17 93 116 1 0 0.83333333 0.03806935 0.8978906 0.77341763 ## 18 95 115 1 0 0.82608696 0.03905833 0.8918099 0.76520757 ## 19 105 114 1 0 0.81884058 0.04003974 0.8856890 0.75703763 ## 20 107 113 1 0 0.81159420 0.04101457 0.8795299 0.74890594 ## 省略一部分 ## 192 550 15 1 0 0.34325958 0.18475887 0.4930486 0.23897674 ## 193 551 14 0 1 0.34325958 0.18475887 0.4930486 0.23897674 ## 194 559 13 0 1 0.34325958 0.18475887 0.4930486 0.23897674 ## 195 588 12 0 1 0.34325958 0.18475887 0.4930486 0.23897674 ## 196 641 11 1 0 0.31205416 0.20791044 0.4690333 0.20761384 ## 197 654 10 1 0 0.28084875 0.23310483 0.4434980 0.17784977 ## 198 687 9 1 0 0.24964333 0.26120251 0.4165393 0.14961805 ## 199 705 8 1 0 0.21843791 0.29340057 0.3882139 0.12290937 ## 200 728 7 1 0 0.18723250 0.33150176 0.3585552 0.09777018 ## 201 731 6 1 0 0.15602708 0.37845310 0.3275969 0.07431220 ## 202 735 5 1 0 0.12482167 0.43957565 0.2954319 0.05273786 ## 203 740 4 0 1 0.12482167 0.43957565 0.2954319 0.05273786 ## 204 765 3 1 0 0.08321444 0.59991117 0.2696771 0.02567754 ## 205 821 2 0 1 0.08321444 0.59991117 0.2696771 0.02567754 ## 206 965 1 0 1 0.08321444 0.59991117 0.2696771 0.02567754 ## strata sex ## 1 sex=1 1 ## 2 sex=1 1 ## 3 sex=1 1 ## 4 sex=1 1 ## 5 sex=1 1 ## 6 sex=1 1 ## 7 sex=1 1 ## 8 sex=1 1 ## 9 sex=1 1 ## 10 sex=1 1 ## 11 sex=1 1 ## 12 sex=1 1 ## 13 sex=1 1 ## 14 sex=1 1 ## 15 sex=1 1 ## 16 sex=1 1 ## 17 sex=1 1 ## 18 sex=1 1 ## 19 sex=1 1 ## 20 sex=1 1 ## 省略一部分 ## 192 sex=2 2 ## 193 sex=2 2 ## 194 sex=2 2 ## 195 sex=2 2 ## 196 sex=2 2 ## 197 sex=2 2 ## 198 sex=2 2 ## 199 sex=2 2 ## 200 sex=2 2 ## 201 sex=2 2 ## 202 sex=2 2 ## 203 sex=2 2 ## 204 sex=2 2 ## 205 sex=2 2 ## 206 sex=2 2 通过ggsurvplot()进行可视化，非常多的细节可以修改，超级详细的教程可以参考后面的推文：超级 详细的R语言生存分析可视化
ggsurvplot(fit.logrank, "hv" , # Add medians survival # Change legends: title & labels "Sex" ,"Male" , "Female" ),# Add p-value and tervals TRUE , # 这里P值直接写数字也行 TRUE ,# Add risk table TRUE , 0.2 ,TRUE ,# Color palettes. Use custom color: c("#E7B800", "#2E9FDF"), # or brewer color (e.g.: "Dark2"), or ggsci color (e.g.: "jco") "#E7B800" , "#2E9FDF" ),# Change ggplot2 theme # Change font size, style and color "Survival curve" ,16 , "bold" , "darkblue" ),14 , "bold.italic" , "red" ),14 , "bold.italic" , "darkred" ),12 , "plain" , "darkgreen" )自带的surv_cutpoint()可用于寻找最佳切点，但是只能用于连续性数据。
使用myeloma数据进行演示。
rm(list = ls())# 0. Load some data ## molecular_group chr1q21_status treatment event time CCND1 CRIM1 ## GSM50986 Cyclin D-1 3 copies TT2 0 69.24 9908.4 420.9 ## GSM50988 Cyclin D-2 2 copies TT2 0 66.43 16698.8 52.0 ## GSM50989 MMSET 2 copies TT2 0 66.50 294.5 617.9 ## GSM50990 MMSET 3 copies TT2 1 42.67 241.9 11.9 ## GSM50991 MAF <NA> TT2 0 65.00 472.6 38.8 ## GSM50992 Hyperdiploid 2 copies TT2 0 65.20 664.1 16.9 ## DEPDC1 IRF4 TP53 WHSC1 ## GSM50986 523.5 16156.5 10.0 261.9 ## GSM50988 21.1 16946.2 1056.9 363.8 ## GSM50989 192.9 8903.9 1762.8 10042.9 ## GSM50990 184.7 11894.7 946.8 4931.0 ## GSM50991 212.0 7563.1 361.4 165.0 ## GSM50992 341.6 16023.4 2096.3 569.2 寻找最佳切点：
# 1. Determine the optimal cutpoint of variables "time" , event = "event" ,"DEPDC1" , "WHSC1" , "CRIM1" ) # 找这3个变量的最佳切点 ## cutpoint statistic ## DEPDC1 279.8 4.275452 ## WHSC1 3205.6 3.361330 ## CRIM1 82.3 1.968317 查看根据最佳切点进行分组后的数据分布情况：
# 2. Plot cutpoint for DEPDC1 "DEPDC1" , palette = "npg" )## $DEPDC1 根据最佳切点重新划分数据，这样数据就根据最佳切点变成了高表达/低表达组。
# 3. Categorize variables ## time event DEPDC1 WHSC1 CRIM1 ## GSM50986 69.24 0 high low high ## GSM50988 66.43 0 low low low ## GSM50989 66.50 0 low high high ## GSM50990 42.67 1 low high low ## GSM50991 65.00 0 low low low ## GSM50992 65.20 0 high low low 根据最佳切点绘制生存曲线：
# 4. Fit survival curves and visualize library ("survival" )TRUE , conf.int = TRUE )确定最佳切点的R包还有非常多，其他的后续会再介绍。
多时间点和多指标的ROC曲线绘制，可参考：R语言画多时间点ROC和多指标ROC曲线
平滑版ORC曲线和最佳截点：生存资料ROC曲线的最佳截点和平滑曲线
ROC曲线的比较：ROC（AUC）曲线的显著性检验
下次继续介绍Cox回归。
参考资料 http://www./english/wiki/survival-analysis-basics https://www./tutorials/survival_analysis_in_r_tutorial.html https://mp.weixin.qq.com/s/2rwxeaF_M0UnqPi2F9JNxA 
