生存分析
梁雪枫
生存分析(Survival analysis)是指根据试验或调查得到的数据对生物或人的生存时间进行分析和推断，研究生存时间和结局与众多影响因素间关系及其程度大小的方法，也称生存率分析或存活率分析。生存分析适合于处理时间－事件数据,生存时间(survival time)是指从某起点事件开始到被观测对象出现终点事件所经历的时间，如从疾病的“确诊”到“死亡”。生存时间有两种类型：完全数据(complete data)指被观测对象从观察起点到出现终点事件所经历的时间;截尾数据(consored data)或删失数据，指在出现终点事件前，被观测对象的观测过程终止了。由于被观测对象所提供的信息是不完全的，只知道他们的生存事件超过了截尾时间。截尾主要由于失访、退出和终止产生。生存分析方法大体上可分为三类：非参数法、半参数方法和参数法，用Kaplan-Meier曲线(也称乘积极限法Product limit method)和寿命表法(Life table method)估计生存率和中位生存时间等是非参数的方法，半参数方法指Cox比例风险模型，参数方法指指数模型、Weibull模型、Gompertz模型等分析方法。
死亡概率(mortality probability)指某段时间开始时生存的个体在该段时间内死亡的可能性大小，若无删失数据，死亡概率=某人群某段时间总死亡例数/该人群同时间段期初观察例数。生存概率(survival probability)指某段时间开始时存活的个体至该时间结束时仍然存活的可能性大小，生存概率=1-死亡概率=某人群活过某段时间例数/该人群同时间段期初观察例数。由于生存分析中常存在删失数据，假定删失事件在观察时间内各个时间点等机会发生，分母改用校正观察例数。校正观察例数=期初观察例数-删失例数/2。 生存率(Survival rate)，用\(S(t_{k})\)表示，指经历\(t_{k}\)个单位时间后仍存活的概率，若无删失数据，则为活过了\(t_{k}\)时刻仍然存活的例数/观察开始的总例数。如果有删失数据，分母则需要按时段进行校正，此时生存率的计算公式为\[S(t_{k})=P(T>t_{k})=p_{1} \cdot p_{2}\cdot \cdot \cdot p_{k}\],其中\(p_{1} \cdot p_{2}\cdot \cdot \cdot p_{k}\)表示不同时间段的生存概率。生存率为多个时间段生存概率的累积，故又称累积生存概率，其标准误计算公式为\[SE(S(t_{k}))=S(t_{k})\sqrt{\sum_{i=1}^{k}\frac{q_{i}}{p_{i}n_{i}}}\],\(q_{i}\)为死亡概率，\(p_{i}\)为生存概率。
例 addicts是238名病例随访信息，status变量表示病例的生存状况（0为删失，1为终点事件），Days.survival变量表示生存的天数。
addicts <- read.table('ADDICTS.txt',T)
addicts$Clinic <- as.factor(addicts$Clinic)
addicts$Prison <- as.factor(addicts$Prison)非参数法
寿命表（Life Table）
寿命表时描述一段时间内生存状况、终点事件和生存概率的表格，需计算累积生存概率即每一步生存概率的乘积，可完成对病例随访资料在任意指定时点的生存状况评价。survival包中包括了所有生存分析所必须的函数，生存分析主要是把数据放入Surv object，通过Surv()函数做进一步分析。Surv object是将时间和生存状况的信息合并在一个简单的对象内，Surv(time, time2, event,type=c(‘right’, ‘left’, ‘interval’, ‘counting’, ‘interval2’, ‘mstate’),origin=0)，time为生存时间，time2为区间删失的结束时间，event为生存状况，生存状况变量必须是数值或者逻辑型的。如果时数值型，则有两个选项，0表示删失，1表示终点事件，或者1表示删失，2表示终点事件。如果时逻辑型的，则FALSE表示删失，True表示终点事件。type为删失的类型有右删失、左删失、区间删失、第一类区间删失、第二类区间删失。
addicts$surv <- Surv(addicts$Days.survival,addicts$Status)
summary(survfit(addicts$surv~1),censor=T)## Call: survfit(formula = addicts$surv ~ 1)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 2 238 0 1.000 0.00000 1.0000 1.000
## 7 236 1 0.996 0.00423 0.9875 1.000
## 13 235 1 0.992 0.00597 0.9799 1.000
## 17 234 1 0.987 0.00729 0.9731 1.000
## 19 233 1 0.983 0.00840 0.9667 1.000
## 26 232 1 0.979 0.00937 0.9606 0.997
## 28 231 0 0.979 0.00937 0.9606 0.997
## 29 229 1 0.975 0.01026 0.9546 0.995
## 30 228 1 0.970 0.01107 0.9488 0.992
## 33 227 1 0.966 0.01182 0.9431 0.989
## 35 226 2 0.957 0.01317 0.9320 0.984
## 37 224 1 0.953 0.01379 0.9265 0.981
## 41 223 2 0.945 0.01493 0.9158 0.974
## 47 221 1 0.940 0.01546 0.9105 0.971
## 49 220 1 0.936 0.01597 0.9053 0.968
## 50 219 1 0.932 0.01646 0.9001 0.965
## 53 218 0 0.932 0.01646 0.9001 0.965
## 59 216 1 0.927 0.01694 0.8949 0.961
## 62 215 1 0.923 0.01740 0.8897 0.958
## 67 213 1 0.919 0.01785 0.8845 0.954
## 72 212 0 0.919 0.01785 0.8845 0.954
## 75 211 1 0.914 0.01829 0.8793 0.951
## 79 210 1 0.910 0.01871 0.8742 0.948
## 84 209 1 0.906 0.01913 0.8691 0.944
## 86 208 0 0.906 0.01913 0.8691 0.944
## 90 207 1 0.901 0.01953 0.8639 0.940
## 95 206 1 0.897 0.01992 0.8588 0.937
## 96 205 1 0.893 0.02029 0.8537 0.933
## 98 204 0 0.893 0.02029 0.8537 0.933
## 103 203 0 0.893 0.02029 0.8537 0.933
## 109 202 1 0.888 0.02067 0.8486 0.930
## 111 201 0 0.888 0.02067 0.8486 0.930
## 117 200 1 0.884 0.02104 0.8435 0.926
## 122 199 1 0.879 0.02140 0.8384 0.922
## 126 198 1 0.875 0.02174 0.8333 0.919
## 127 197 1 0.870 0.02208 0.8282 0.915
## 129 196 1 0.866 0.02241 0.8232 0.911
## 136 194 1 0.862 0.02274 0.8181 0.907
## 143 193 1 0.857 0.02305 0.8131 0.903
## 145 192 1 0.853 0.02336 0.8080 0.900
## 146 191 0 0.853 0.02336 0.8080 0.900
## 147 190 1 0.848 0.02366 0.8030 0.896
## 148 189 0 0.848 0.02366 0.8030 0.896
## 149 188 1 0.844 0.02396 0.7979 0.892
## 150 187 1 0.839 0.02426 0.7929 0.888
## 157 185 1 0.835 0.02455 0.7878 0.884
## 160 184 1 0.830 0.02483 0.7828 0.880
## 161 183 1 0.826 0.02510 0.7777 0.876
## 167 181 1 0.821 0.02538 0.7727 0.872
## 168 180 1 0.816 0.02564 0.7676 0.868
## 170 179 1 0.812 0.02590 0.7626 0.864
## 175 178 1 0.807 0.02615 0.7576 0.860
## 176 176 1 0.803 0.02640 0.7526 0.856
## 180 175 2 0.794 0.02689 0.7425 0.848
## 181 173 1 0.789 0.02712 0.7375 0.844
## 183 172 1 0.784 0.02735 0.7325 0.840
## 190 171 1 0.780 0.02757 0.7275 0.836
## 192 170 1 0.775 0.02779 0.7226 0.832
## 193 169 1 0.771 0.02800 0.7176 0.827
## 204 168 1 0.766 0.02821 0.7127 0.823
## 205 166 1 0.761 0.02841 0.7077 0.819
## 207 165 1 0.757 0.02861 0.7027 0.815
## 209 164 1 0.752 0.02881 0.6978 0.811
## 210 163 0 0.752 0.02881 0.6978 0.811
## 212 162 2 0.743 0.02919 0.6878 0.802
## 216 160 2 0.734 0.02955 0.6779 0.794
## 222 158 0 0.734 0.02955 0.6779 0.794
## 223 157 1 0.729 0.02973 0.6729 0.790
## 231 156 1 0.724 0.02991 0.6679 0.785
## 232 155 1 0.720 0.03008 0.6630 0.781
## 237 154 1 0.715 0.03024 0.6580 0.777
## 244 153 1 0.710 0.03040 0.6531 0.772
## 247 152 1 0.706 0.03056 0.6481 0.768
## 257 151 1 0.701 0.03071 0.6432 0.764
## 258 150 1 0.696 0.03086 0.6383 0.759
## 259 149 1 0.692 0.03101 0.6333 0.755
## 262 148 2 0.682 0.03128 0.6235 0.746
## 268 146 2 0.673 0.03154 0.6138 0.738
## 275 144 1 0.668 0.03167 0.6089 0.733
## 280 143 1 0.663 0.03179 0.6040 0.729
## 283 142 0 0.663 0.03179 0.6040 0.729
## 286 141 1 0.659 0.03191 0.5991 0.724
## 293 140 1 0.654 0.03203 0.5942 0.720
## 294 139 1 0.649 0.03214 0.5893 0.716
## 299 138 1 0.645 0.03225 0.5844 0.711
## 302 137 1 0.640 0.03236 0.5796 0.707
## 314 136 1 0.635 0.03246 0.5747 0.702
## 317 135 0 0.635 0.03246 0.5747 0.702
## 322 134 1 0.631 0.03256 0.5698 0.698
## 325 133 0 0.631 0.03256 0.5698 0.698
## 326 132 0 0.631 0.03256 0.5698 0.698
## 337 131 1 0.626 0.03267 0.5648 0.693
## 341 129 1 0.621 0.03277 0.5598 0.689
## 342 128 0 0.621 0.03277 0.5598 0.689
## 346 127 0 0.621 0.03277 0.5598 0.689
## 348 126 1 0.616 0.03288 0.5547 0.684
## 350 125 1 0.611 0.03298 0.5496 0.679
## 358 124 1 0.606 0.03308 0.5446 0.675
## 366 122 1 0.601 0.03318 0.5395 0.670
## 367 121 1 0.596 0.03328 0.5343 0.665
## 368 119 1 0.591 0.03338 0.5292 0.660
## 376 118 1 0.586 0.03347 0.5241 0.656
## 386 117 1 0.581 0.03355 0.5189 0.651
## 389 116 1 0.576 0.03364 0.5138 0.646
## 393 115 1 0.571 0.03371 0.5087 0.641
## 394 114 1 0.566 0.03379 0.5036 0.636
## 399 112 1 0.561 0.03386 0.4984 0.631
## 405 111 0 0.561 0.03386 0.4984 0.631
## 408 110 0 0.561 0.03386 0.4984 0.631
## 428 109 1 0.556 0.03394 0.4932 0.627
## 434 108 1 0.551 0.03401 0.4879 0.622
## 438 107 1 0.546 0.03408 0.4827 0.617
## 439 106 0 0.546 0.03408 0.4827 0.617
## 450 105 1 0.540 0.03415 0.4774 0.612
## 452 104 1 0.535 0.03422 0.4722 0.607
## 456 103 0 0.535 0.03422 0.4722 0.607
## 457 102 1 0.530 0.03428 0.4668 0.602
## 460 101 1 0.525 0.03434 0.4615 0.597
## 461 100 0 0.525 0.03434 0.4615 0.597
## 465 99 1 0.519 0.03440 0.4562 0.591
## 475 98 0 0.519 0.03440 0.4562 0.591
## 480 97 0 0.519 0.03440 0.4562 0.591
## 482 96 1 0.514 0.03447 0.4507 0.586
## 489 95 1 0.509 0.03453 0.4452 0.581
## 496 94 1 0.503 0.03458 0.4398 0.576
## 504 92 1 0.498 0.03463 0.4342 0.570
## 512 91 1 0.492 0.03468 0.4287 0.565
## 514 90 1 0.487 0.03473 0.4232 0.560
## 517 89 1 0.481 0.03476 0.4178 0.554
## 518 87 1 0.476 0.03480 0.4122 0.549
## 522 86 1 0.470 0.03483 0.4067 0.544
## 523 85 2 0.459 0.03488 0.3956 0.533
## 531 83 0 0.459 0.03488 0.3956 0.533
## 532 80 1 0.453 0.03491 0.3899 0.527
## 533 78 1 0.448 0.03495 0.3841 0.522
## 540 77 1 0.442 0.03497 0.3783 0.516
## 541 76 0 0.442 0.03497 0.3783 0.516
## 543 75 0 0.442 0.03497 0.3783 0.516
## 546 74 1 0.436 0.03501 0.3723 0.510
## 550 73 1 0.430 0.03503 0.3664 0.504
## 551 72 0 0.430 0.03503 0.3664 0.504
## 555 71 0 0.430 0.03503 0.3664 0.504
## 560 70 1 0.424 0.03507 0.3603 0.498
## 563 69 1 0.418 0.03509 0.3542 0.492
## 564 66 0 0.418 0.03509 0.3542 0.492
## 566 64 0 0.418 0.03509 0.3542 0.492
## 575 63 0 0.418 0.03509 0.3542 0.492
## 581 62 1 0.411 0.03517 0.3474 0.486
## 587 60 0 0.411 0.03517 0.3474 0.486
## 591 59 1 0.404 0.03525 0.3404 0.479
## 602 57 0 0.404 0.03525 0.3404 0.479
## 609 56 0 0.404 0.03525 0.3404 0.479
## 611 55 0 0.404 0.03525 0.3404 0.479
## 612 54 2 0.389 0.03550 0.3252 0.465
## 613 52 0 0.389 0.03550 0.3252 0.465
## 624 51 1 0.381 0.03561 0.3175 0.458
## 633 50 0 0.381 0.03561 0.3175 0.458
## 641 49 0 0.381 0.03561 0.3175 0.458
## 646 48 1 0.373 0.03574 0.3095 0.450
## 652 47 1 0.365 0.03586 0.3015 0.443
## 661 46 1 0.357 0.03595 0.2935 0.435
## 667 45 1 0.350 0.03601 0.2856 0.428
## 679 44 1 0.342 0.03606 0.2777 0.420
## 683 43 1 0.334 0.03609 0.2699 0.412
## 684 41 0 0.334 0.03609 0.2699 0.412
## 708 39 1 0.325 0.03616 0.2614 0.404
## 713 38 0 0.325 0.03616 0.2614 0.404
## 714 37 1 0.316 0.03623 0.2527 0.396
## 730 36 0 0.316 0.03623 0.2527 0.396
## 739 35 1 0.307 0.03631 0.2437 0.387
## 749 34 1 0.298 0.03635 0.2348 0.379
## 755 33 1 0.289 0.03635 0.2260 0.370
## 760 32 1 0.280 0.03632 0.2173 0.361
## 769 31 0 0.280 0.03632 0.2173 0.361
## 771 28 1 0.270 0.03638 0.2075 0.352
## 774 27 1 0.260 0.03638 0.1978 0.342
## 785 26 1 0.250 0.03633 0.1882 0.332
## 787 25 0 0.250 0.03633 0.1882 0.332
## 788 24 0 0.250 0.03633 0.1882 0.332
## 790 23 0 0.250 0.03633 0.1882 0.332
## 796 22 0 0.250 0.03633 0.1882 0.332
## 808 21 0 0.250 0.03633 0.1882 0.332
## 821 20 2 0.225 0.03675 0.1635 0.310
## 826 18 0 0.225 0.03675 0.1635 0.310
## 836 17 1 0.212 0.03690 0.1506 0.298
## 837 16 1 0.199 0.03689 0.1380 0.286
## 840 15 0 0.199 0.03689 0.1380 0.286
## 857 14 1 0.184 0.03688 0.1246 0.273
## 878 13 1 0.170 0.03667 0.1116 0.260
## 881 12 0 0.170 0.03667 0.1116 0.260
## 884 11 0 0.170 0.03667 0.1116 0.260
## 892 10 1 0.153 0.03675 0.0958 0.245
## 899 9 1 0.136 0.03639 0.0807 0.230
## 905 8 0 0.136 0.03639 0.0807 0.230
## 932 7 0 0.136 0.03639 0.0807 0.230
## 944 5 0 0.136 0.03639 0.0807 0.230
## 969 4 0 0.136 0.03639 0.0807 0.230
## 1021 3 0 0.136 0.03639 0.0807 0.230
## 1052 2 0 0.136 0.03639 0.0807 0.230
## 1076 1 0 0.136 0.03639 0.0807 0.230上表的第一行表示，在第2天，有238个调查对象，没有发生终点事件（n.event），生存概率为(238-0)/238=1,其中有2个删失对象没有显示出来。第二行表示在第7天，有236个对象，其中1个发生了终点事件，生存概率为(236-1)/235*1=0.996。寿命表中其他数据行的意思类似。
Kaplan-Meier曲线
Kaplan-Meier曲线也称生存曲线， 纵轴表示生存概率，横轴表示生存事件，它是一条下降的曲线,下降的坡度越陡,表示生存率越低或生存时间越短,其斜率表示死亡速率。如果在概率50%处画一条水平线，它将中位生存事件点和生存曲线相交。
KM0 <- survfit(surv ~ 1, type="kaplan-meier",data=addicts)
kml <- summary(KM0,censor=T)
attributes(kml)## $names
## [1] "n" "time" "n.risk" "n.event" "n.censor"
## [6] "surv" "type" "std.err" "upper" "lower"
## [11] "conf.type" "conf.int" "call" "table"
##
## $class
## [1] "summary.survfit"plot(kml$time,kml$surv,type="s")
plot(survfit(addicts$surv~1))
#绘制一条曲线时，图形中a包含95%的置信区间和删失标记，不需要可设置为False
plot(survfit(addicts$surv~1),conf.int = F,mark.time = F)
abline(h=0.5,lty=2,col="red")
#中位生存期
survfit(addicts$surv~1)## Call: survfit(formula = addicts$surv ~ 1)
##
## n events median 0.95LCL 0.95UCL
## 238 150 504 399 560#25%,50%和75%生存期
quantile(KM0, probs=c(0.25, 0.5, 0.75), conf.int=FALSE)## 25 50 75
## 212 504 821#50天和100天生存状况
summary(KM0, times=c(50, 100))## Call: survfit(formula = surv ~ 1, data = addicts, type = "kaplan-meier")
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 50 219 16 0.932 0.0165 0.900 0.965
## 100 203 9 0.893 0.0203 0.854 0.933在survfit函数中改变公式右边的参数，可获得不同因子水平的生存曲线。
plot(survfit(addicts$surv~addicts$Clinic),col=c("red","blue"),conf.int = F)
legend(10,.4,legend=c("1","2"),col = c("red","blue"),lty=c(1,1))
不同生存曲线间是否有差异，可通过survdiff进行比较，该函数最后一个参数时rho，用于指定检验的类型。让rho=0(默认时)，进行对数秩(log-rank)检验或\(Mantel-Haenszel \chi^{2}\)检验，比较各组期望频数和实际观察数。如果两组间的差异水平太大，\(\chi^{2}\)会较大而P值较小，表示生存曲线有统计学差异。当rho=1时，进行Gehan-Wilcoxon的Peto校正检验，该检验赋予早期终点事件较大的权重。
survdiff(addicts$surv~addicts$Clinic)## Call:
## survdiff(formula = addicts$surv ~ addicts$Clinic)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## addicts$Clinic=1 163 122 90.9 10.6 27.9
## addicts$Clinic=2 75 28 59.1 16.4 27.9
##
## Chisq= 27.9 on 1 degrees of freedom, p= 1.28e-07分层比较
在Clinic变量有可能和其他变量之间存在相关性，应调整其影响后，研究生存区间之间的差异。
cc(addicts$Clinic,addicts$Prison)
##
## addicts$Prison
## addicts$Clinic 0 1 Total
## 1 88 75 163
## 2 39 36 75
## Total 127 111 238
##
## OR = 1.08
## Exact 95% CI = 0.6, 1.94
## Chi-squared = 0.08, 1 d.f., P value = 0.775
## Fisher's exact test (2-sided) P value = 0.782survdiff(addicts$surv~addicts$Clinic+strata(addicts$Prison))## Call:
## survdiff(formula = addicts$surv ~ addicts$Clinic + strata(addicts$Prison))
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## addicts$Clinic=1 163 122 91.7 10.0 26.9
## addicts$Clinic=2 75 28 58.3 15.8 26.9
##
## Chisq= 26.9 on 1 degrees of freedom, p= 2.1e-07在调整addicts$Prison的影响之后，与原始情况没有太大差异，说明Prison对surv的影响不是独立的。
累积风险率
风险率指每个单位时间的时小比例，这个随时间变化。用图形可绘制累积风险率，其斜率可相对容易的观察。
plot(survfit(addicts$surv~1),conf.int = F,fun="cumhaz")
上图显示，在后800多天的时，由于没有终点事件的发生，斜率时水平的。
在survfit函数中改变公式右边的参数，可获得不同因子水平的累积风险概率。
plot(survfit(addicts$surv~addicts$Clinic),col=c("red","blue"),conf.int = F,fun="cumhaz")
legend(5,4,legend=c("1","2"),col = c("red","blue"),lty=c(1,1))
model <-coxph(surv~Clinic,data=addicts)包含通过寿命表（Life Table）分析法，；Kaplan-Meier方法，对病例随访资料进行生存分析，在对应于每一实际观察事件时点上，作生存率的评价和建立Cox回归模型（亦称比例风险模型）。
参数法(Parametric proportional hazards models)
参数方法要求观察的生存事件服从某一特定的分布，采用估计分布中参数的方法获得生存率的估计值。生存事件的分布可能为指数分布、weibull分布、对数正态分布等，这些分布曲线都有相应的生存率函数形式，只需求的相应参数的估计值，即可获得生存率的估计值和生存曲线。
假定生存时间符合weibull分布
fitWeib <- survreg(surv~Clinic+Prison+Dose, dist="weibull", data=addicts)
summary(fitWeib)##
## Call:
## survreg(formula = surv ~ Clinic + Prison + Dose, data = addicts,
## dist = "weibull")
## Value Std. Error z p
## (Intercept) 4.8139 0.27499 17.51 1.29e-68
## Clinic2 0.7090 0.15722 4.51 6.49e-06
## Prison1 -0.2295 0.12079 -1.90 5.75e-02
## Dose 0.0244 0.00459 5.32 1.03e-07
## Log(scale) -0.3150 0.06756 -4.66 3.13e-06
##
## Scale= 0.73
##
## Weibull distribution
## Loglik(model)= -1084.5 Loglik(intercept only)= -1114.9
## Chisq= 60.89 on 3 degrees of freedom, p= 3.8e-13
## Number of Newton-Raphson Iterations: 7
## n= 238AFT参数转换为Cox模型的\(\beta\)
(betaHat <- -coef(fitWeib) / fitWeib$scale)## (Intercept) Clinic2 Prison1 Dose
## -6.59596008 -0.97152449 0.31441429 -0.03346751模型比较
fitExp <- survreg(surv~Clinic+Prison+Dose,dist="exponential", data=addicts)
anova(fitExp, fitWeib) ## Terms Resid. Df -2*LL Test Df Deviance Pr(>Chi)
## 1 Clinic + Prison + Dose 234 2187.942 NA NA NA
## 2 Clinic + Prison + Dose 233 2168.953 = 1 18.98925 1.314567e-05提出因子变量后，模型的比较
fitR <- survreg(surv~Dose, dist="weibull", data=addicts)
anova(fitR, fitWeib)## Terms Resid. Df -2*LL Test Df Deviance Pr(>Chi)
## 1 Dose 235 2195.360 NA NA NA
## 2 Clinic + Prison + Dose 233 2168.953 = 2 26.40651 1.844589e-06生存曲线估计
dfNew <- data.frame(Clinic=factor(c("1", "2"), levels=levels(addicts$Clinic)),
Dose=c(50, 60),
Prison=factor(c("0", "1"), levels=levels(addicts$Prison)))
percs <- (1:99)/100
FWeib <- predict(fitWeib, newdata=dfNew, type="quantile", p=percs, se=TRUE)
matplot(cbind(FWeib$fit[1, ],
FWeib$fit[1, ] - 2*FWeib$se.fit[1, ],
FWeib$fit[1, ] + 2*FWeib$se.fit[1, ]), 1-percs,
type="l", main=expression(paste("Weibull-Fit ", hat(S)(t), " mit SE")),
xlab="t", ylab="Survival", lty=c(1, 2, 2), lwd=2, col="blue")
matlines(cbind(FWeib$fit[2, ],
FWeib$fit[2, ] - 2*FWeib$se.fit[2, ],
FWeib$fit[2, ] + 2*FWeib$se.fit[2, ]), 1-percs, col="red", lwd=2)
legend(x="topright", lwd=2, lty=c(1, 2, 1, 2), col=c("blue", "blue", "red", "red"),
legend=c("Clinic=1, Dose=50, Prison=0", "+- 2*SE", "Clinic=2, Dose=60, Prison=1", "+- 2*SE"))
半参数法(COX回归)
多数生存时间的分布并不符合指数分布、weibull分布等，不宜采用参数法进行分析。COX回归用于研究各种因素（称为协变量）对于生存期长短的关系，Cox 回归是一种半参数模型，只规定了影响因素和生存时间的关系，但是没有对生存时间的分布情况加以限定，与参数模型相比，该模型不能给出各时点的风险率，，但可估计出各研究因素对风险率的影响，进行多因素分析。 风险函数(Hazard Function),用h(t)表示,其定义为:\[h(t)=\lim_{\Delta t\rightarrow 0}\left [ \frac{S(t)-S(t+\Delta t)}{\Delta t} \right ]/S(t)\],表示时刻t上一个事件瞬时发生的概率，即一个到t时刻存活的个体，在t时刻事件的瞬时发生率。Cox模型为 \[\ln h(t) = \ln h_{0}(t) + \beta_{1} X_{1} + \dots + \beta_{p} X_{p}\] 其中\(X_{1} + \dots + X_{p}\)是协变量，\(\beta_{1} + \dots + \beta_{p}\)是回归系数，由样本估计而得。\(\beta_{i}>0\)表示该协变量是危险因素，越大使生存时间越短,\(\beta_{i}<0\)表示该协变量是保护因素，越大使生存时间越长。\(h_{0}(t)\)为基础风险函数，它是全部协变量都为0或标准状态下的风险函数，一般是未知的。\(h(t)\)表示当各协变量值X固定时的风险函数，它和\(h_{0}(t)\)成比例，所以该模型又称为比例风险模型（proportional hazard model),COX回归模型不用于估计生存率，主要用于因素分析。
model <- coxph(surv~Clinic+Prison,data=addicts)
summary(model)## Call:
## coxph(formula = surv ~ Clinic + Prison, data = addicts)
##
## n= 238, number of events= 150
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Clinic2 -1.1089 0.3299 0.2144 -5.173 2.3e-07 ***
## Prison1 0.2778 1.3203 0.1656 1.678 0.0933 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Clinic2 0.3299 3.0309 0.2168 0.5022
## Prison1 1.3203 0.7574 0.9544 1.8263
##
## Concordance= 0.599 (se = 0.025 )
## Rsquare= 0.132 (max possible= 0.997 )
## Likelihood ratio test= 33.78 on 2 df, p=4.625e-08
## Wald test = 28.17 on 2 df, p=7.643e-07
## Score (logrank) test = 30.51 on 2 df, p=2.366e-07Clinic2系数为负值，且有统计学意义。Clinic2exp(coef)是0.3412，提示Clinic2与Clinic1相比风险率降低了65.88%(1-0.3412)。Prison1无统计学意义。
模型拟合
AIC值
extractAIC(model)## [1] 2.0 1381.3McFadden, Cox & Snell and Nagelkerke pseudo \(R^{2}\)
LLf <- model$loglik[2]
LL0 <- model$loglik[1]McFadden pseudo-\(R^2\)
as.vector(1 - (LLf / LL0))## [1] 0.02393795Cox & Snell
as.vector(1 - exp((2/nrow(addicts)) * (LL0 - LLf)))## [1] 0.1323143Nagelkerke
as.vector((1 - exp((2/nrow(addicts)) * (LL0 - LLf))) / (1 - exp(LL0)^(2/nrow(addicts))))## [1] 0.1326674模型比较
model1 <- coxph(surv~Clinic,data=addicts)
anova(model1, model) ## Analysis of Deviance Table
## Cox model: response is surv
## Model 1: ~ Clinic
## Model 2: ~ Clinic + Prison
## loglik Chisq Df P(>|Chi|)
## 1 -690.04
## 2 -688.65 2.7887 1 0.09493 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1生存函数(Survival function)
(CPH <- survfit(model))## Call: survfit(formula = model)
##
## n events median 0.95LCL 0.95UCL
## 238 150 518 450 591quantile(CPH, probs=c(0.25, 0.5, 0.75), conf.int=FALSE)## 25 50 75
## 231 518 821估计生存函数(Estimated survival function for new specific data)
dfNew <- data.frame(Clinic=factor(c("2", "2"), levels=levels(addicts$Clinic)),
X=c(-2, -2),
Prison=factor(c("0", "1"), levels=levels(addicts$Prison)))
CPHnew <- survfit(model, newdata=dfNew)
par(mar=c(5, 4.5, 4, 2)+0.1, cex.lab=1.4, cex.main=1.4)
plot(CPH, main=expression(paste("Cox PH-estimate ", hat(S)(t), " with CI")),
xlab="t", ylab="Survival", lwd=2)
lines(CPHnew$time, CPHnew$surv[ , 1], lwd=2, col="blue")
lines(CPHnew$time, CPHnew$surv[ , 2], lwd=2, col="red")
legend(x="topright", lwd=2, col=c("black", "blue", "red"),
legend=c("pseudo-observation", "Clinic=2, X=-2, Prison=0", "Clinic=2, X=-2, Prison=1"))
累积基础风险函数(Cumulative baseline hazard)
expCoef <- exp(coef(model))
Lambda0A <- basehaz(model, centered=FALSE)
Lambda0B <- expCoef[2]*Lambda0A$hazard
Lambda0C <- expCoef[3]*Lambda0A$hazard
plot(hazard ~ time, main=expression(paste("Cox PH-estimate ", hat(Lambda)[g](t), " per group")),
type="s", ylim=c(0, 5), xlab="t", ylab="cumulative hazard", lwd=2, data=Lambda0A)
lines(Lambda0A$time, Lambda0B, lwd=2, col="red")
lines(Lambda0A$time, Lambda0C, lwd=2, col="green")
legend(x="bottomright", lwd=2, col=1:3, legend=LETTERS[1:3])
模型诊断(Model diagnostics)
比例风险假定(Proportional hazards assumption)
将纵轴取对数后，绘制时间的对数值图形，可以比较Clinic变量两种取值的生存曲线。如果两条曲线平行，则不太可能违反比例风险假定。
plot(survfit(surv~Clinic,data=addicts),fun="cloglog",conf.int = F,col = c("red","blue"))
两条曲线相交不止一次，从图形很难判断是否违反比例风险假定。可采取如下检验
czph <- cox.zph(coxph(surv~Clinic+Prison,data=addicts))
czph ## rho chisq p
## Clinic2 -0.2698 12.179 0.000483
## Prison1 -0.0292 0.128 0.720595
## GLOBAL NA 12.663 0.001779par(mfrow=c(2, 2))
plot(czph)
结果显示，违反比例风险假定的证据非常强。图形展现的是随时间变化的\(\beta\)图形。
影响分析(Influential observations)
dfbetas <- residuals(coxph(surv~Clinic+Prison,data=addicts), type="dfbetas")
par(mfrow=c(1, 2))
plot(dfbetas[ , 1], type="h", main="DfBETAS for Clinic", ylab="DfBETAS", lwd=2)
plot(dfbetas[ , 2], type="h", main="DfBETAS for Prison", ylab="DfBETAS", lwd=2)
对数线性假设(Linearity of log hazard)
对自变量是连续性变量，需检测其线性的假设
resMart <- residuals(coxph(surv~Clinic+Dose,data=addicts), type="martingale")
par(mfrow=c(1, 1))
plot(addicts$Dose, resMart, main="Martingale-residuals for Dose",
xlab="Dose", ylab="Residuen", pch=20)
lines(loess.smooth(addicts$Dose, resMart), lwd=2, col="blue")
legend(x="bottomleft", col="blue", lwd=2, legend="LOESS fit")
预测风险(Predicted hazard ratios)
根据建立的模型对每个个体的风险率进行预测，连续性变量假定其等于样本均值，因子变量假定其等于亚变量。 ####风险率预测
predRes <- predict(coxph(surv~Clinic+Prison,data=addicts), type="risk")
head(predRes, n=10) #显示前10个## [1] 1.245897 1.644901 1.245897 1.245897 1.644901 1.245897 1.644901
## [8] 1.644901 1.245897 1.644901生存期预测
Shat1 <- survexp(~ 1, ratetable=model, data=addicts)
with(Shat1, head(data.frame(time, surv), n=4))## time surv
## 1 7 0.9957744
## 2 13 0.9915437
## 3 17 0.9873255
## 4 19 0.9831022分因子变量生存期预测
Shat2 <- survexp(~ Clinic, ratetable=model, data=addicts)
with(Shat2, head(data.frame(time, surv), n=4))## time Clinic.1 Clinic.2
## 1 7 0.9946485 0.9982213
## 2 13 0.9892931 0.9964348
## 3 17 0.9839563 0.9946479
## 4 19 0.9786156 0.9928531分层回归
phFit <- coxph(surv~Clinic+Prison+Dose,data=addicts)
summary(phFit)## Call:
## coxph(formula = surv ~ Clinic + Prison + Dose, data = addicts)
##
## n= 238, number of events= 150
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Clinic2 -1.009896 0.364257 0.214889 -4.700 2.61e-06 ***
## Prison1 0.326555 1.386184 0.167225 1.953 0.0508 .
## Dose -0.035369 0.965249 0.006379 -5.545 2.94e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Clinic2 0.3643 2.7453 0.2391 0.5550
## Prison1 1.3862 0.7214 0.9988 1.9238
## Dose 0.9652 1.0360 0.9533 0.9774
##
## Concordance= 0.665 (se = 0.026 )
## Rsquare= 0.238 (max possible= 0.997 )
## Likelihood ratio test= 64.56 on 3 df, p=6.228e-14
## Wald test = 54.12 on 3 df, p=1.056e-11
## Score (logrank) test = 56.32 on 3 df, p=3.598e-12step(phFit)## Start: AIC=1352.52
## surv ~ Clinic + Prison + Dose
##
## Df AIC
## <none> 1352.5
## - Prison 1 1354.3
## - Clinic 1 1376.9
## - Dose 1 1381.3## Call:
## coxph(formula = surv ~ Clinic + Prison + Dose, data = addicts)
##
##
## coef exp(coef) se(coef) z p
## Clinic2 -1.00990 0.36426 0.21489 -4.70 2.6e-06
## Prison1 0.32655 1.38618 0.16722 1.95 0.051
## Dose -0.03537 0.96525 0.00638 -5.54 2.9e-08
##
## Likelihood ratio test=64.6 on 3 df, p=6.23e-14
## n= 238, number of events= 150没有变量剔除时，AIC水平是最低的，所有变量均应保留。
cox.zph(phFit)## rho chisq p
## Clinic2 -0.2578 11.19 0.000824
## Prison1 -0.0382 0.22 0.639369
## Dose 0.0724 0.70 0.402749
## GLOBAL NA 12.62 0.005546全局检验的P值，有统计学意义，表面违反了比例风险的假定。一种可能方法是对因子变量Clinic进行分层分析。
strataphFit <- coxph(surv~strata(Clinic)+Prison+Dose,data=addicts)
cox.zph(strataphFit)## rho chisq p
## Prison1 -0.0205 0.0628 0.802
## Dose 0.0860 0.9953 0.318
## GLOBAL NA 1.0186 0.601summary(strataphFit)## Call:
## coxph(formula = surv ~ strata(Clinic) + Prison + Dose, data = addicts)
##
## n= 238, number of events= 150
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Prison1 0.389605 1.476397 0.168930 2.306 0.0211 *
## Dose -0.035115 0.965495 0.006465 -5.432 5.59e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Prison1 1.4764 0.6773 1.0603 2.0559
## Dose 0.9655 1.0357 0.9533 0.9778
##
## Concordance= 0.651 (se = 0.034 )
## Rsquare= 0.133 (max possible= 0.994 )
## Likelihood ratio test= 33.91 on 2 df, p=4.322e-08
## Wald test = 32.66 on 2 df, p=8.076e-08
## Score (logrank) test = 33.33 on 2 df, p=5.774e-08用Clinic进行分层降低了\(\chi ^{2}\)值，全局检验的P值，没有统计学意义，没有违反比列风险的假定。strataphFit模型与phFit模型比较，Clinic分层因素的系统被忽略。