|
(文章中公式可以查看博文http://www./post/r-survival/)
學(xué)習(xí)生存分析預(yù)先要求對(duì)R有所了解,，基本能夠操作R數(shù)據(jù)框和包的使用,。要是懂ggplot2和dplyr就更好了,。
資料
-
生存分析備查表 - 會(huì)使用到的主要函數(shù)與包介紹
-
背景介紹 - 統(tǒng)計(jì)知識(shí)
背景
在縱向研究中,，我們需要從一個(gè)時(shí)間點(diǎn)追蹤樣本或受試者（例如,，進(jìn)入研究，診斷,，開始治療）,，直到我們觀察到一些結(jié)果事件（例如死亡，疾病發(fā)作，復(fù)發(fā)）,，但不會(huì)有意義的假設(shè)改變的速率是不變的,。例如：手術(shù)后心臟手術(shù)后的死亡風(fēng)險(xiǎn)最高，術(shù)后恢復(fù)的患者死亡風(fēng)險(xiǎn)緩慢降低,，隨著患者年齡的增長(zhǎng),，風(fēng)險(xiǎn)再次緩慢上升?；蛘?，不同癌癥的復(fù)發(fā)率隨時(shí)間變化很大，并且取決于腫瘤遺傳學(xué),，治療和其他環(huán)境因素,。
定義
生存分析可讓我們分析事件發(fā)生的速率，而不會(huì)假設(shè)速率不變,。一般而言,，生存分析可以讓我們對(duì)事件發(fā)生之前的時(shí)間進(jìn)行建模1或比較不同組之間的事件時(shí)間，或者事件時(shí)間與定量變量之間的相關(guān)性,。
風(fēng)險(xiǎn)是特定時(shí)間點(diǎn)t的瞬時(shí)事件發(fā)生（死亡）率,。生存分析并不認(rèn)為隨著時(shí)間的推移危害是不變的。累積風(fēng)險(xiǎn)是直到時(shí)間t為止經(jīng)歷的總風(fēng)險(xiǎn),。生存函數(shù)是個(gè)體在時(shí)間t之前存在的概率（或者,，不發(fā)生感興趣事件的概率）。這是事件（例如,，死亡）尚未發(fā)生的可能性,。看起來像這樣,，其中TT是死亡時(shí)間,，Pr(T>t)Pr(T>t)是死亡時(shí)間大于某個(gè)時(shí)間tt的概率,。SS是概率,，所以0≤S(t)≤10≤S(t)≤1，因?yàn)樯嫫诳偸钦担═≥0T≥0）,。
S(t)=Pr(T>t)S(t)=Pr(T>t)
Kaplan-Meier曲線描述了生存函數(shù),。這是一個(gè)階梯函數(shù)，說明隨著時(shí)間的推移累計(jì)生存概率,。曲線在沒有事件發(fā)生的時(shí)間段內(nèi)是水平的,，然后垂直下降，對(duì)應(yīng)于每次發(fā)生事件時(shí)生存函數(shù)的變化,。截尾是一種生存分析特有的缺失數(shù)據(jù)問題,。 當(dāng)我們?cè)谘芯拷Y(jié)束時(shí)跟蹤樣本/主題并且事件從未發(fā)生時(shí)會(huì)發(fā)生這種情況。這也可能是由于樣本/受試者因死亡以外的原因而退出研究或其他一些失訪導(dǎo)致的。樣本數(shù)據(jù)發(fā)生截尾是因?yàn)槟阒恢肋@個(gè)人存活到失去跟蹤為止,，但你不知道任何關(guān)于之后他的生存狀態(tài)2,。
比例風(fēng)險(xiǎn)假設(shè)：生存分析的主要目標(biāo)是比較不同組別（例如白血病患者與正常對(duì)照組）的生存功能。如果兩組人都死亡,，那么兩條生存曲線都將結(jié)束于0％,，但是其中一組的平均存活時(shí)間可能比另一組長(zhǎng)。生存分析通過比較觀察期間不同時(shí)間的風(fēng)險(xiǎn)來做到這一點(diǎn),。生存分析并不假設(shè)危害是恒定的,，但確實(shí)假定組間危害的比率隨著時(shí)間的推移是恒定的。3本文不包括處理非比例風(fēng)險(xiǎn)的方法或伴隨時(shí)間到事件的協(xié)變量交互作用,。
比例風(fēng)險(xiǎn)回歸也稱為Cox回歸,，是評(píng)估不同變量對(duì)生存率影響的最常用方法。
Cox PH模型
Kaplan-Meier曲線適用于觀察兩個(gè)分類組4之間生存率的差異,，但對(duì)于評(píng)估諸如年齡,，基因表達(dá)，白細(xì)胞計(jì)數(shù)等定量變量的影響,，它們不起作用,。Cox PH回歸可評(píng)估分類變量和連續(xù)變量的影響，并且可以一次模擬多個(gè)變量的影響,。
Cox PH回歸模型將時(shí)間t處的風(fēng)險(xiǎn)自然對(duì)數(shù)表示為h(t)h(t),，作為基線危險(xiǎn)（h0(t)h0(t)）的函數(shù)（所有暴露變量為0的個(gè)體的風(fēng)險(xiǎn)）和多個(gè)暴露變量x1x1，x1x1,，......,，xpxp。 Cox PH模型的形式是：
log(h(t))=log(h0(t))+β1x1+β2x2+...+βpxplog(h(t))=log(h0(t))+β1x1+β2x2+...+βpxp
如果對(duì)方程的兩邊進(jìn)行冪運(yùn)算,，并將右邊限制為僅包含兩個(gè)組（x1=1x1=1作為暴露變量,，x1=0x1=0作為非暴露變量）的單個(gè)分類暴露變量（x1x1），則等式變?yōu)椋?/p>
h1(t)=h0(t)×eβ1x1h1(t)=h0(t)×eβ1x1
重新排列該等式可以估計(jì)風(fēng)險(xiǎn)比率,，比較時(shí)間t暴露對(duì)于未暴露個(gè)體：
HR(t)=h1(t)h0(t)=eβ1HR(t)=h1(t)h0(t)=eβ1
該模型顯示風(fēng)險(xiǎn)比為eβ1eβ1,，并且在時(shí)間t內(nèi)保持不變（因此名稱比例風(fēng)險(xiǎn)回歸）。ββ值是根據(jù)模型估計(jì)的回歸系數(shù),，并表示相應(yīng)預(yù)測(cè)變量中每單位增加的log(Hazard,Ratio)log(Hazard,Ratio),。危害比的解釋取決于預(yù)測(cè)變量的測(cè)量尺度，但簡(jiǎn)單地說,，正系數(shù)表示較差的存活率,，負(fù)系數(shù)表示所討論的變量存活率較高。
用R進(jìn)行生存分析
核心的分析函數(shù)都在survival包里,，我們這里使用dplyr包,，然后用survminer包進(jìn)行繪圖,。
# 確保在導(dǎo)入前安裝好
library(survival)
library(dplyr)
library(survminer)
我們將使用的核心函數(shù)包括：
-
Surv()：創(chuàng)建一個(gè)生存對(duì)象
-
survfit()：使用公式或已構(gòu)建的Cox模型擬合生存曲線
-
coxph()：擬合Cox比例風(fēng)險(xiǎn)回歸模型
其他我們可能會(huì)用到的函數(shù)：
-
cox.zph()：檢驗(yàn)一個(gè)Cox回歸模型的比例風(fēng)險(xiǎn)假設(shè)
-
survdiff()：用log-rank/Mantel-Haenszel檢驗(yàn)檢驗(yàn)生存差異5
Surv()創(chuàng)建響應(yīng)變量，典型用法是使用事件時(shí)間,，6以及事件是否發(fā)生（即死亡與截尾）,。survfit()創(chuàng)建一條生存曲線，然后可以顯示或繪圖,。coxph()實(shí)現(xiàn)回歸分析,，并且模型以與常規(guī)線性模型中相同的方式指定，但使用coxph()函數(shù),。
開始
我們將使用內(nèi)置的肺癌數(shù)據(jù)集7學(xué)習(xí)使用survival包,。你可以通過運(yùn)行?lung獲取數(shù)據(jù)集信息：
library(survival)
?lung
-
inst: Institution code
-
time: Survival time in days
-
status: censoring status 1=censored, 2=dead
-
age: Age in years
-
sex: Male=1 Female=2
-
ph.ecog: ECOG performance score (0=good 5=dead)
-
ph.karno: Karnofsky performance score as rated by physician
-
pat.karno: Karnofsky performance score as rated by patient
-
meal.cal: Calories consumed at meals
-
wt.loss: Weight loss in last six months
lung <- as_tibble(lung)
lung
## # A tibble: 228 x 10
## inst time status age sex ph.ecog ph.karno pat.karno meal.cal
## * <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3. 306. 2. 74. 1. 1. 90. 100. 1175.
## 2 3. 455. 2. 68. 1. 0. 90. 90. 1225.
## 3 3. 1010. 1. 56. 1. 0. 90. 90. NA
## 4 5. 210. 2. 57. 1. 1. 90. 60. 1150.
## 5 1. 883. 2. 60. 1. 0. 100. 90. NA
## 6 12. 1022. 1. 74. 1. 1. 50. 80. 513.
## 7 7. 310. 2. 68. 2. 2. 70. 60. 384.
## 8 11. 361. 2. 71. 2. 2. 60. 80. 538.
## 9 1. 218. 2. 53. 1. 1. 70. 80. 825.
## 10 7. 166. 2. 61. 1. 2. 70. 70. 271.
## # ... with 218 more rows, and 1 more variable: wt.loss <dbl>
生存曲線
構(gòu)建生存對(duì)象：
s <- Surv(time = lung$time, event = lung$status)
class(s)
## [1] "Surv"
s
## [1] 306 455 1010+ 210 883 1022+ 310 361 218 166 170
## [12] 654 728 71 567 144 613 707 61 88 301 81
## [23] 624 371 394 520 574 118 390 12 473 26 533
## [34] 107 53 122 814 965+ 93 731 460 153 433 145
## [45] 583 95 303 519 643 765 735 189 53 246 689
## [56] 65 5 132 687 345 444 223 175 60 163 65
## [67] 208 821+ 428 230 840+ 305 11 132 226 426 705
## [78] 363 11 176 791 95 196+ 167 806+ 284 641 147
## [89] 740+ 163 655 239 88 245 588+ 30 179 310 477
## [100] 166 559+ 450 364 107 177 156 529+ 11 429 351
## [111] 15 181 283 201 524 13 212 524 288 363 442
## [122] 199 550 54 558 207 92 60 551+ 543+ 293 202
## [133] 353 511+ 267 511+ 371 387 457 337 201 404+ 222
## [144] 62 458+ 356+ 353 163 31 340 229 444+ 315+ 182
## [155] 156 329 364+ 291 179 376+ 384+ 268 292+ 142 413+
## [166] 266+ 194 320 181 285 301+ 348 197 382+ 303+ 296+
## [177] 180 186 145 269+ 300+ 284+ 350 272+ 292+ 332+ 285
## [188] 259+ 110 286 270 81 131 225+ 269 225+ 243+ 279+
## [199] 276+ 135
## [ reached getOption("max.print") -- omitted 28 entries ]
現(xiàn)在，讓我們使用survfit()函數(shù)擬合一條生存曲線,。這里讓我們先創(chuàng)建一條不考慮任何比較的生存曲線,，所以我們只需要指定survfit()在公式里期望的截距（比如~1。
survfit(s~1)
## Call: survfit(formula = s ~ 1)
##
## n events median 0.95LCL 0.95UCL
## 228 165 310 285 363
# 前面操作可以一步完成
survfit(Surv(time, status)~1, data=lung)
## Call: survfit(formula = Surv(time, status) ~ 1, data = lung)
##
## n events median 0.95LCL 0.95UCL
## 228 165 310 285 363
但模型對(duì)象本身不會(huì)給出太多的價(jià)值信息,，我們需要使用summary函數(shù)查看模型匯總結(jié)果：
sfit <- survfit(Surv(time, status)~1, data=lung)
summary(sfit)
## Call: survfit(formula = Surv(time, status) ~ 1, data = lung)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 5 228 1 0.9956 0.00438 0.9871 1.000
## 11 227 3 0.9825 0.00869 0.9656 1.000
## 12 224 1 0.9781 0.00970 0.9592 0.997
## 13 223 2 0.9693 0.01142 0.9472 0.992
## 15 221 1 0.9649 0.01219 0.9413 0.989
## 26 220 1 0.9605 0.01290 0.9356 0.986
## 30 219 1 0.9561 0.01356 0.9299 0.983
## 31 218 1 0.9518 0.01419 0.9243 0.980
## 53 217 2 0.9430 0.01536 0.9134 0.974
## 54 215 1 0.9386 0.01590 0.9079 0.970
## 59 214 1 0.9342 0.01642 0.9026 0.967
## 60 213 2 0.9254 0.01740 0.8920 0.960
## 61 211 1 0.9211 0.01786 0.8867 0.957
## 62 210 1 0.9167 0.01830 0.8815 0.953
## 65 209 2 0.9079 0.01915 0.8711 0.946
## 71 207 1 0.9035 0.01955 0.8660 0.943
## 79 206 1 0.8991 0.01995 0.8609 0.939
## 81 205 2 0.8904 0.02069 0.8507 0.932
## 88 203 2 0.8816 0.02140 0.8406 0.925
## 92 201 1 0.8772 0.02174 0.8356 0.921
## 93 199 1 0.8728 0.02207 0.8306 0.917
## 95 198 2 0.8640 0.02271 0.8206 0.910
## 105 196 1 0.8596 0.02302 0.8156 0.906
## 107 194 2 0.8507 0.02362 0.8056 0.898
## 110 192 1 0.8463 0.02391 0.8007 0.894
## 116 191 1 0.8418 0.02419 0.7957 0.891
## 118 190 1 0.8374 0.02446 0.7908 0.887
## 122 189 1 0.8330 0.02473 0.7859 0.883
## [到達(dá)getOption("max.print") -- 略過111行]]
這個(gè)表格每一行顯示了一個(gè)（多個(gè)）事件或截尾發(fā)生了,，在風(fēng)險(xiǎn)中的樣本數(shù)（就是還沒死的），以及及時(shí)的累積生存率等,。
如果我們不使用截距建模,，結(jié)果會(huì)更加有意思：
sfit <- survfit(Surv(time, status)~sex, data=lung)
sfit
## Call: survfit(formula = Surv(time, status) ~ sex, data = lung)
##
## n events median 0.95LCL 0.95UCL
## sex=1 138 112 270 212 310
## sex=2 90 53 426 348 550
summary(sfit)
## Call: survfit(formula = Surv(time, status) ~ sex, data = lung)
##
## sex=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 11 138 3 0.9783 0.0124 0.9542 1.000
## 12 135 1 0.9710 0.0143 0.9434 0.999
## 13 134 2 0.9565 0.0174 0.9231 0.991
## 15 132 1 0.9493 0.0187 0.9134 0.987
## 26 131 1 0.9420 0.0199 0.9038 0.982
## 30 130 1 0.9348 0.0210 0.8945 0.977
## 31 129 1 0.9275 0.0221 0.8853 0.972
## 53 128 2 0.9130 0.0240 0.8672 0.961
## 54 126 1 0.9058 0.0249 0.8583 0.956
## 59 125 1 0.8986 0.0257 0.8496 0.950
## 60 124 1 0.8913 0.0265 0.8409 0.945
## 65 123 2 0.8768 0.0280 0.8237 0.933
## 71 121 1 0.8696 0.0287 0.8152 0.928
## 81 120 1 0.8623 0.0293 0.8067 0.922
## 88 119 2 0.8478 0.0306 0.7900 0.910
## 92 117 1 0.8406 0.0312 0.7817 0.904
## 93 116 1 0.8333 0.0317 0.7734 0.898
## 95 115 1 0.8261 0.0323 0.7652 0.892
## 105 114 1 0.8188 0.0328 0.7570 0.886
## 107 113 1 0.8116 0.0333 0.7489 0.880
## 110 112 1 0.8043 0.0338 0.7408 0.873
## 116 111 1 0.7971 0.0342 0.7328 0.867
## 118 110 1 0.7899 0.0347 0.7247 0.861
## 131 109 1 0.7826 0.0351 0.7167 0.855
## 132 108 2 0.7681 0.0359 0.7008 0.842
## 135 106 1 0.7609 0.0363 0.6929 0.835
## 142 105 1 0.7536 0.0367 0.6851 0.829
## 144 104 1 0.7464 0.0370 0.6772 0.823
## [到達(dá)getOption("max.print") -- 略過71行]]
##
## sex=2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 5 90 1 0.9889 0.0110 0.9675 1.000
## 60 89 1 0.9778 0.0155 0.9478 1.000
## 61 88 1 0.9667 0.0189 0.9303 1.000
## 62 87 1 0.9556 0.0217 0.9139 0.999
## 79 86 1 0.9444 0.0241 0.8983 0.993
## 81 85 1 0.9333 0.0263 0.8832 0.986
## 95 83 1 0.9221 0.0283 0.8683 0.979
## 107 81 1 0.9107 0.0301 0.8535 0.972
## 122 80 1 0.8993 0.0318 0.8390 0.964
## 145 79 2 0.8766 0.0349 0.8108 0.948
## 153 77 1 0.8652 0.0362 0.7970 0.939
## 166 76 1 0.8538 0.0375 0.7834 0.931
## 167 75 1 0.8424 0.0387 0.7699 0.922
## 182 71 1 0.8305 0.0399 0.7559 0.913
## 186 70 1 0.8187 0.0411 0.7420 0.903
## 194 68 1 0.8066 0.0422 0.7280 0.894
## 199 67 1 0.7946 0.0432 0.7142 0.884
## 201 66 2 0.7705 0.0452 0.6869 0.864
## 208 62 1 0.7581 0.0461 0.6729 0.854
## 226 59 1 0.7452 0.0471 0.6584 0.843
## 239 57 1 0.7322 0.0480 0.6438 0.833
## 245 54 1 0.7186 0.0490 0.6287 0.821
## 268 51 1 0.7045 0.0501 0.6129 0.810
## 285 47 1 0.6895 0.0512 0.5962 0.798
## 293 45 1 0.6742 0.0523 0.5791 0.785
## 305 43 1 0.6585 0.0534 0.5618 0.772
## 310 42 1 0.6428 0.0544 0.5447 0.759
## 340 39 1 0.6264 0.0554 0.5267 0.745
## [到達(dá)getOption("max.print") -- 略過23行]]
summary()函數(shù)中可以設(shè)定時(shí)間參數(shù)用來選定一個(gè)時(shí)間區(qū)間，我們可以以此比對(duì)男生是不是比女生有更高的風(fēng)險(xiǎn)：
summary(sfit, times=seq(0, 1000, 100))
## Call: survfit(formula = Surv(time, status) ~ sex, data = lung)
##
## sex=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 138 0 1.0000 0.0000 1.0000 1.000
## 100 114 24 0.8261 0.0323 0.7652 0.892
## 200 78 30 0.6073 0.0417 0.5309 0.695
## 300 49 20 0.4411 0.0439 0.3629 0.536
## 400 31 15 0.2977 0.0425 0.2250 0.394
## 500 20 7 0.2232 0.0402 0.1569 0.318
## 600 13 7 0.1451 0.0353 0.0900 0.234
## 700 8 5 0.0893 0.0293 0.0470 0.170
## 800 6 2 0.0670 0.0259 0.0314 0.143
## 900 2 2 0.0357 0.0216 0.0109 0.117
## 1000 2 0 0.0357 0.0216 0.0109 0.117
##
## sex=2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 90 0 1.0000 0.0000 1.0000 1.000
## 100 82 7 0.9221 0.0283 0.8683 0.979
## 200 66 11 0.7946 0.0432 0.7142 0.884
## 300 43 9 0.6742 0.0523 0.5791 0.785
## 400 26 10 0.5089 0.0603 0.4035 0.642
## 500 21 5 0.4110 0.0626 0.3050 0.554
## 600 11 3 0.3433 0.0634 0.2390 0.493
## 700 8 3 0.2496 0.0652 0.1496 0.417
## 800 2 5 0.0832 0.0499 0.0257 0.270
## 900 1 0 0.0832 0.0499 0.0257 0.270
Kaplan-Meier生存曲線
現(xiàn)在我們使用Kaplan-Meier曲線來可視化這一結(jié)果：
plot(sfit)
R的plot()函數(shù)選項(xiàng)可以用來修改這個(gè)圖,，你可以參加?plot.survfit,。我們這里不會(huì)描述太多細(xì)節(jié)，因?yàn)橛辛硪粋€(gè)叫survminer的包提供的一個(gè)叫ggsurvplot()的函數(shù)可以幫助我們更簡(jiǎn)單地做出可以發(fā)表的生存曲線,，如果你對(duì)ggplot2語法很熟悉的話還能更簡(jiǎn)單地進(jìn)行修改,。讓我們導(dǎo)入并嘗試一下吧：
library(survminer)
ggsurvplot(sfit)
這個(gè)圖比剛才那個(gè)圖更好看，信息量也更多：它用顏色幫我們區(qū)分了組別,，并添加了橫縱坐標(biāo)的標(biāo)簽,。
讓我們添加曲線的置信區(qū)間，并增加long-rank檢驗(yàn)的結(jié)果p值以及風(fēng)險(xiǎn)表格：
ggsurvplot(sfit, conf.int=TRUE, pval=TRUE, risk.table=TRUE,
legend.labs=c("Male", "Female"), legend.title="Sex",
palette=c("dodgerblue2", "orchid2"),
title="Kaplan-Meier Curve for Lung Cancer Survival",
risk.table.height=.15)
Cox回歸模型
Kaplan-Meier曲線用來對(duì)兩個(gè)分類變量差異的可視化非常合適,，但分類要是多,，那就糟透了：
ggsurvplot(survfit(Surv(time, status)~nodes, data=survival::colon))
而且生存曲線另外不能可視化的是連續(xù)型變量的風(fēng)險(xiǎn)。
Cox PH回歸模型正好是處理這類問題的一把好手,，它同樣內(nèi)置于survival包中,，語法與lm()和glm()一致,。
讓我們?cè)賮碛梅伟?shù)據(jù)集看看不同性別的風(fēng)險(xiǎn),，這次使用Cox模型。
fit <- coxph(Surv(time, status)~sex, data=lung)
fit
## Call:
## coxph(formula = Surv(time, status) ~ sex, data = lung)
##
## coef exp(coef) se(coef) z p
## sex -0.531 0.588 0.167 -3.18 0.0015
##
## Likelihood ratio test=10.6 on 1 df, p=0.00111
## n= 228, number of events= 165
結(jié)果中的exp(coef)列包含eβ1eβ1,。它就是風(fēng)險(xiǎn)比率——該變量對(duì)風(fēng)險(xiǎn)率的乘數(shù)效應(yīng)（對(duì)于該變量每個(gè)單位增加的）,。因此,，對(duì)于像性別這樣的分類變量，從男性（基線）到女性的結(jié)果大約減少約40％的危險(xiǎn),。你也可以翻轉(zhuǎn)coef列上的符號(hào),，并采用exp(0.531)，你可以將其解釋為男性導(dǎo)致危險(xiǎn)增加1.7倍,，或者單位時(shí)間男性的死亡率約為女性1.7倍（女性死亡率為男性的0.588倍）,。
要記住：
- HR=1: 沒有效應(yīng)
- HR>1: 風(fēng)險(xiǎn)增加
- HR<1: 風(fēng)險(xiǎn)減少 （保護(hù)變量）
你還會(huì)注意到，“性別”有一個(gè)對(duì)應(yīng)的p值,，整個(gè)模型中也有一個(gè)p值,。0.00111的p值非常接近我們?cè)贙aplan-Meier圖上看到的p=0.00131的p值。這是因?yàn)镵M曲線顯示的是對(duì)數(shù)秩檢驗(yàn)的p值,。你可以通過調(diào)用summary(fit)來獲得Cox模型結(jié)果,。你也可以使用survdiff()直接計(jì)算log-rank測(cè)試p值。
summary(fit)
## Call:
## coxph(formula = Surv(time, status) ~ sex, data = lung)
##
## n= 228, number of events= 165
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sex -0.531 0.588 0.167 -3.18 0.0015 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex 0.588 1.7 0.424 0.816
##
## Concordance= 0.579 (se = 0.022 )
## Rsquare= 0.046 (max possible= 0.999 )
## Likelihood ratio test= 10.6 on 1 df, p=0.00111
## Wald test = 10.1 on 1 df, p=0.00149
## Score (logrank) test = 10.3 on 1 df, p=0.00131
survdiff(Surv(time, status)~sex, data=lung)
## Call:
## survdiff(formula = Surv(time, status) ~ sex, data = lung)
##
## 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.00131
讓我們創(chuàng)建另一個(gè)模型,，分析數(shù)據(jù)集中的所有變量,！這向我們展示了當(dāng)所有變量一起考慮時(shí)，如何影響生存,。一些是非常強(qiáng)大的預(yù)測(cè)指標(biāo)（性別,，ECOG評(píng)分）。有趣的是,，醫(yī)師對(duì)Karnofsky表現(xiàn)評(píng)分的評(píng)分稍高,，但患者評(píng)分相同。
fit <- coxph(Surv(time, status)~sex+age+ph.ecog+ph.karno+pat.karno+meal.cal+wt.loss, data=lung)
fit
## Call:
## coxph(formula = Surv(time, status) ~ sex + age + ph.ecog + ph.karno +
## pat.karno + meal.cal + wt.loss, data = lung)
##
## coef exp(coef) se(coef) z p
## sex -5.51e-01 5.76e-01 2.01e-01 -2.74 0.0061
## age 1.06e-02 1.01e+00 1.16e-02 0.92 0.3591
## ph.ecog 7.34e-01 2.08e+00 2.23e-01 3.29 0.0010
## ph.karno 2.25e-02 1.02e+00 1.12e-02 2.00 0.0457
## pat.karno -1.24e-02 9.88e-01 8.05e-03 -1.54 0.1232
## meal.cal 3.33e-05 1.00e+00 2.60e-04 0.13 0.8979
## wt.loss -1.43e-02 9.86e-01 7.77e-03 -1.84 0.0652
##
## Likelihood ratio test=28.3 on 7 df, p=0.000192
## n= 168, number of events= 121
## (60 observations deleted due to missingness)
分類畫KM曲線
讓我們回到肺部數(shù)據(jù)并查看年齡的Cox模型,?？雌饋砟挲g在模擬為連續(xù)變量時(shí)似乎有一點(diǎn)重要。
coxph(Surv(time, status)~age, data=lung)
## Call:
## coxph(formula = Surv(time, status) ~ age, data = lung)
##
## coef exp(coef) se(coef) z p
## age 0.0187 1.0189 0.0092 2.03 0.042
##
## Likelihood ratio test=4.24 on 1 df, p=0.0395
## n= 228, number of events= 165
現(xiàn)在我們的的回歸分析顯示年齡有重要意義,，讓我們制作Kaplan-Meier圖,。但是，正如我們之前所看到的,，我們不能這樣做,，因?yàn)槲覀儠?huì)為每個(gè)獨(dú)特的年齡值獲得單獨(dú)的曲線！
ggsurvplot(survfit(Surv(time, status)~age, data=lung))
你可能在這里看到的一件事是試圖將一個(gè)連續(xù)變量分成不同的組 - 三分位數(shù),，上四分位數(shù)與下四分位數(shù),，中位數(shù)分?jǐn)?shù)等 - 這樣你就可以生成生存曲線圖。但是,，你如何進(jìn)行分組是有意義的,！檢查cut的幫助。cut()接受一個(gè)連續(xù)變量和一些斷點(diǎn),，并從中創(chuàng)建一個(gè)分類變量,。 我們來得到數(shù)據(jù)集的平均年齡,，并繪制一個(gè)顯示年齡分布的直方圖。
mean(lung$age)
hist(lung$age)
ggplot(lung, aes(age)) + geom_histogram(bins=20)
現(xiàn)在,，讓我們嘗試通過lung$age創(chuàng)建一個(gè)分類變量,，其中0,62（平均值）和正無窮大。我們可以在這里繼續(xù)添加labels =選項(xiàng)來標(biāo)記我們創(chuàng)建的分組,，例如,，“年輕”和“老”。最后,，我們可以將這個(gè)結(jié)果分配給肺數(shù)據(jù)集中的一個(gè)新對(duì)象,。
cut(lung$age, breaks=c(0, 62, Inf))
## [1] (62,Inf] (62,Inf] (0,62] (0,62] (0,62] (62,Inf] (62,Inf]
## [8] (62,Inf] (0,62] (0,62] (0,62] (62,Inf] (62,Inf] (0,62]
## [15] (0,62] (62,Inf] (62,Inf] (62,Inf] (0,62] (0,62] (62,Inf]
## [22] (0,62] (0,62] (0,62] (62,Inf] (62,Inf] (0,62] (62,Inf]
## [29] (0,62] (62,Inf] (62,Inf] (62,Inf] (0,62] (0,62] (0,62]
## [36] (0,62] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (62,Inf]
## [43] (0,62] (0,62] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (62,Inf]
## [50] (0,62] (62,Inf] (62,Inf] (62,Inf] (0,62] (0,62] (0,62]
## [57] (62,Inf] (0,62] (0,62] (62,Inf] (62,Inf] (0,62] (62,Inf]
## [64] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (62,Inf]
## [71] (62,Inf] (0,62] (62,Inf] (0,62] (0,62] (62,Inf] (0,62]
## [78] (0,62] (62,Inf] (62,Inf] (0,62] (0,62] (0,62] (0,62]
## [85] (0,62] (62,Inf] (0,62] (0,62] (0,62] (62,Inf] (62,Inf]
## [92] (62,Inf] (62,Inf] (0,62] (62,Inf] (62,Inf] (62,Inf] (62,Inf]
## [99] (62,Inf] (62,Inf] (0,62] (62,Inf] (0,62] (62,Inf] (0,62]
## [106] (62,Inf] (0,62] (62,Inf] (0,62] (62,Inf] (62,Inf] (0,62]
## [113] (62,Inf] (62,Inf] (0,62] (62,Inf] (0,62] (62,Inf] (62,Inf]
## [120] (62,Inf] (62,Inf] (0,62] (62,Inf] (62,Inf] (62,Inf] (62,Inf]
## [127] (0,62] (62,Inf] (62,Inf] (0,62] (0,62] (0,62] (0,62]
## [134] (0,62] (62,Inf] (62,Inf] (0,62] (0,62] (0,62] (0,62]
## [141] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (0,62] (0,62] (62,Inf]
## [148] (0,62] (62,Inf] (0,62] (62,Inf] (0,62] (0,62] (0,62]
## [155] (0,62] (62,Inf] (62,Inf] (0,62] (62,Inf] (0,62] (0,62]
## [162] (0,62] (62,Inf] (62,Inf] (62,Inf] (0,62] (0,62] (0,62]
## [169] (0,62] (62,Inf] (0,62] (0,62] (0,62] (0,62] (0,62]
## [176] (0,62] (0,62] (0,62] (0,62] (62,Inf] (0,62] (0,62]
## [183] (62,Inf] (62,Inf] (0,62] (0,62] (62,Inf] (0,62] (62,Inf]
## [190] (0,62] (62,Inf] (0,62] (0,62] (62,Inf] (62,Inf] (62,Inf]
## [197] (62,Inf] (62,Inf] (0,62] (0,62]
## [ reached getOption("max.print") -- omitted 28 entries ]
## Levels: (0,62] (62,Inf]
cut(lung$age, breaks=c(0, 62, Inf), labels=c("young", "old"))
## [1] old old young young young old old old young young young
## [12] old old young young old old old young young old young
## [23] young young old old young old young old old old young
## [34] young young young old old old old old old young young
## [45] old old old old old young old old old young young
## [56] young old young young old old young old old old old
## [67] old old old old old young old young young old young
## [78] young old old young young young young young old young young
## [89] young old old old old young old old old old old
## [100] old young old young old young old young old young old
## [111] old young old old young old young old old old old
## [122] young old old old old young old old young young young
## [133] young young old old young young young young old old old
## [144] old young young old young old young old young young young
## [155] young old old young old young young young old old old
## [166] young young young young old young young young young young young
## [177] young young young old young young old old young young old
## [188] young old young old young young old old old old old
## [199] young young
## [ reached getOption("max.print") -- omitted 28 entries ]
## Levels: young old
# the base r way:
lung$agecat <- cut(lung$age, breaks=c(0, 62, Inf), labels=c("young", "old"))
# or the dplyr way:
lung <- lung %>%
mutate(agecat=cut(age, breaks=c(0, 62, Inf), labels=c("young", "old")))
head(lung)
## # A tibble: 6 x 11
## inst time status age sex ph.ecog ph.karno pat.karno meal.cal
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3. 306. 2. 74. 1. 1. 90. 100. 1175.
## 2 3. 455. 2. 68. 1. 0. 90. 90. 1225.
## 3 3. 1010. 1. 56. 1. 0. 90. 90. NA
## 4 5. 210. 2. 57. 1. 1. 90. 60. 1150.
## 5 1. 883. 2. 60. 1. 0. 100. 90. NA
## 6 12. 1022. 1. 74. 1. 1. 50. 80. 513.
## # ... with 2 more variables: wt.loss <dbl>, agecat <fct>
現(xiàn)在，當(dāng)我們用這個(gè)新的分類生成KM圖時(shí)會(huì)發(fā)生什么,？ 看起來“老”和“年輕”患者之間的曲線存在一些差異,，老年患者的生存幾率稍差。但是p=0.39時(shí),，62歲以下和62歲以上者的生存率差異不顯著,。
ggsurvplot(survfit(Surv(time, status)~agecat, data=lung), pval=TRUE)
但是，如果我們選擇一個(gè)不同的切點(diǎn),，例如70歲,，這大概是年齡分布上限的四分之一（見“分位數(shù)”）。結(jié)果現(xiàn)在非常重要,！
# the base r way:
lung$agecat <- cut(lung$age, breaks=c(0, 70, Inf), labels=c("young", "old"))
# or the dplyr way:
lung <- lung %>%
mutate(agecat=cut(age, breaks=c(0, 70, Inf), labels=c("young", "old")))
# plot!
ggsurvplot(survfit(Surv(time, status)~agecat, data=lung), pval=TRUE)
請(qǐng)記住,，Cox回歸分析整個(gè)分布范圍內(nèi)的連續(xù)變量，其中Kaplan-Meier圖上的對(duì)數(shù)秩檢驗(yàn)可能會(huì)根據(jù)您對(duì)連續(xù)變量進(jìn)行分類而發(fā)生變化,。他們以一種不同的方式回答類似的問題：回歸模型提出的問題是“年齡對(duì)生存有什么影響,？”，而對(duì)數(shù)秩檢驗(yàn)和KM圖則問：“那些不到70歲和70歲以上的人有差異嗎,？“,。
在新的survminer 0.2.4版本中，新增了可以一次確定一個(gè)或多個(gè)連續(xù)變量最佳分割點(diǎn)的函數(shù)surv_cutpoint()與surv_categorize(),。參閱這篇博文學(xué)習(xí)詳細(xì)的函數(shù)用法,。
英文來源：http:///r-survival.html。上面講述的內(nèi)容涵蓋了R生存分析的基本各個(gè)方面,，具體的使用,、圖形優(yōu)化還需要讀者自己琢磨。英文來源中有不少練習(xí),，有興趣的不妨試試,。后續(xù)還有一些TCGA的分析實(shí)例，有空我再寫寫,。
- 在醫(yī)學(xué)界,，我們通常會(huì)從字面上思考生存分析 - 追蹤直至死亡的時(shí)間,。 但是,，它比這更普遍-生存分析模擬事件發(fā)生之前的時(shí)間（任何事件）都可以,。這可以是生物有機(jī)體的死亡。但也可能是機(jī)械系統(tǒng)發(fā)生硬件故障,，恢復(fù)時(shí)間,，失去工作后有人失業(yè)的時(shí)間，直到成熟的番茄被放牧的鹿吃掉的時(shí)間,，直到有人在車間里睡著的時(shí)間,， 生存分析還包括工程可靠性理論，經(jīng)濟(jì)學(xué)持續(xù)時(shí)間分析和社會(huì)學(xué)事件歷史分析,。?
- 這描述了最常見的截尾類型 - 右截尾,。當(dāng)“開始”未知時(shí)，例如當(dāng)初始診斷或暴露時(shí)間未知時(shí),，通常不會(huì)發(fā)生左側(cè)截尾,。?
- 而且，按照上述定義,，假定兩組之間的累積危險(xiǎn)比率隨著時(shí)間的推移保持不變,。?
- 對(duì)于這些差異，有一個(gè)類似卡方的統(tǒng)計(jì)檢驗(yàn),，稱為對(duì)數(shù)秩檢驗(yàn),，比較生存函數(shù)分類組。?
- Cox回歸和來自
survdiff的logrank檢驗(yàn)將在大多數(shù)時(shí)間給你類似的結(jié)果,。對(duì)數(shù)秩檢驗(yàn)是詢問兩組患者生存曲線是否顯著不同,。Cox回歸是詢問許多分類或連續(xù)變量中哪一個(gè)顯著影響生存。?
-
Surv()也可以輸入開始與截止時(shí)間,，參見?Surv?
- Loprinzi et al. Prospective evaluation of prognostic variables from patient-completed questionnaires. North Central Cancer Treatment Group. Journal of Clinical Oncology. 12(3):601-7, 1994.?
|
