|
(文章中公式可以查看博文http://www./post/r-survival/)
學習生存分析預先要求對R有所了解，基本能夠操作R數(shù)據(jù)框和包的使用。要是懂ggplot2和dplyr就更好了。
資料
背景
在縱向研究中，我們需要從一個時間點追蹤樣本或受試者（例如，進入研究，診斷，開始治療），直到我們觀察到一些結果事件（例如死亡，疾病發(fā)作，復發(fā)），但不會有意義的假設改變的速率是不變的。例如：手術后心臟手術后的死亡風險最高，術后恢復的患者死亡風險緩慢降低，隨著患者年齡的增長，風險再次緩慢上升?；蛘?，不同癌癥的復發(fā)率隨時間變化很大，并且取決于腫瘤遺傳學，治療和其他環(huán)境因素。
定義
生存分析可讓我們分析事件發(fā)生的速率，而不會假設速率不變。一般而言，生存分析可以讓我們對事件發(fā)生之前的時間進行建模1或比較不同組之間的事件時間，或者事件時間與定量變量之間的相關性。
風險是特定時間點t的瞬時事件發(fā)生（死亡）率。生存分析并不認為隨著時間的推移危害是不變的。累積風險是直到時間t為止經歷的總風險。生存函數(shù)是個體在時間t之前存在的概率（或者，不發(fā)生感興趣事件的概率）。這是事件（例如，死亡）尚未發(fā)生的可能性?？雌饋硐襁@樣，其中TT是死亡時間，Pr(T>t)Pr(T>t)是死亡時間大于某個時間tt的概率。SS是概率，所以0≤S(t)≤10≤S(t)≤1，因為生存期總是正值（T≥0T≥0）。
S(t)=Pr(T>t)S(t)=Pr(T>t)
Kaplan-Meier曲線描述了生存函數(shù)。這是一個階梯函數(shù)，說明隨著時間的推移累計生存概率。曲線在沒有事件發(fā)生的時間段內是水平的，然后垂直下降，對應于每次發(fā)生事件時生存函數(shù)的變化。截尾是一種生存分析特有的缺失數(shù)據(jù)問題。 當我們在研究結束時跟蹤樣本/主題并且事件從未發(fā)生時會發(fā)生這種情況。這也可能是由于樣本/受試者因死亡以外的原因而退出研究或其他一些失訪導致的。樣本數(shù)據(jù)發(fā)生截尾是因為你只知道這個人存活到失去跟蹤為止，但你不知道任何關于之后他的生存狀態(tài)2。
比例風險假設：生存分析的主要目標是比較不同組別（例如白血病患者與正常對照組）的生存功能。如果兩組人都死亡，那么兩條生存曲線都將結束于0％，但是其中一組的平均存活時間可能比另一組長。生存分析通過比較觀察期間不同時間的風險來做到這一點。生存分析并不假設危害是恒定的，但確實假定組間危害的比率隨著時間的推移是恒定的。3本文不包括處理非比例風險的方法或伴隨時間到事件的協(xié)變量交互作用。
比例風險回歸也稱為Cox回歸，是評估不同變量對生存率影響的最常用方法。
Cox PH模型
Kaplan-Meier曲線適用于觀察兩個分類組4之間生存率的差異，但對于評估諸如年齡，基因表達，白細胞計數(shù)等定量變量的影響，它們不起作用。Cox PH回歸可評估分類變量和連續(xù)變量的影響，并且可以一次模擬多個變量的影響。
Cox PH回歸模型將時間t處的風險自然對數(shù)表示為h(t)h(t)，作為基線危險（h0(t)h0(t)）的函數(shù)（所有暴露變量為0的個體的風險）和多個暴露變量x1x1，x1x1，......，xpxp。 Cox PH模型的形式是：
log(h(t))=log(h0(t))+β1x1+β2x2+...+βpxplog(h(t))=log(h0(t))+β1x1+β2x2+...+βpxp
如果對方程的兩邊進行冪運算，并將右邊限制為僅包含兩個組（x1=1x1=1作為暴露變量，x1=0x1=0作為非暴露變量）的單個分類暴露變量（x1x1），則等式變?yōu)椋?/p>
h1(t)=h0(t)×eβ1x1h1(t)=h0(t)×eβ1x1
重新排列該等式可以估計風險比率，比較時間t暴露對于未暴露個體：
HR(t)=h1(t)h0(t)=eβ1HR(t)=h1(t)h0(t)=eβ1
該模型顯示風險比為eβ1eβ1，并且在時間t內保持不變（因此名稱比例風險回歸）。ββ值是根據(jù)模型估計的回歸系數(shù)，并表示相應預測變量中每單位增加的log(Hazard,Ratio)log(Hazard,Ratio)。危害比的解釋取決于預測變量的測量尺度，但簡單地說，正系數(shù)表示較差的存活率，負系數(shù)表示所討論的變量存活率較高。
用R進行生存分析
核心的分析函數(shù)都在survival包里，我們這里使用dplyr包，然后用survminer包進行繪圖。
# 確保在導入前安裝好
library(survival)
library(dplyr)
library(survminer)
我們將使用的核心函數(shù)包括：
-
Surv()：創(chuàng)建一個生存對象
-
survfit()：使用公式或已構建的Cox模型擬合生存曲線
-
coxph()：擬合Cox比例風險回歸模型
其他我們可能會用到的函數(shù)：
-
cox.zph()：檢驗一個Cox回歸模型的比例風險假設
-
survdiff()：用log-rank/Mantel-Haenszel檢驗檢驗生存差異5
Surv()創(chuàng)建響應變量，典型用法是使用事件時間，6以及事件是否發(fā)生（即死亡與截尾）。survfit()創(chuàng)建一條生存曲線，然后可以顯示或繪圖。coxph()實現(xiàn)回歸分析，并且模型以與常規(guī)線性模型中相同的方式指定，但使用coxph()函數(shù)。
開始
我們將使用內置的肺癌數(shù)據(jù)集7學習使用survival包。你可以通過運行?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>
生存曲線
構建生存對象：
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
但模型對象本身不會給出太多的價值信息，我們需要使用summary函數(shù)查看模型匯總結果：
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
## [到達getOption("max.print") -- 略過111行]]
這個表格每一行顯示了一個（多個）事件或截尾發(fā)生了，在風險中的樣本數(shù)（就是還沒死的），以及及時的累積生存率等。
如果我們不使用截距建模，結果會更加有意思：
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
## [到達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
## [到達getOption("max.print") -- 略過23行]]
summary()函數(shù)中可以設定時間參數(shù)用來選定一個時間區(qū)間，我們可以以此比對男生是不是比女生有更高的風險：
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曲線來可視化這一結果：
plot(sfit)
R的plot()函數(shù)選項可以用來修改這個圖，你可以參加?plot.survfit。我們這里不會描述太多細節(jié)，因為有另一個叫survminer的包提供的一個叫ggsurvplot()的函數(shù)可以幫助我們更簡單地做出可以發(fā)表的生存曲線，如果你對ggplot2語法很熟悉的話還能更簡單地進行修改。讓我們導入并嘗試一下吧：
library(survminer)
ggsurvplot(sfit)
這個圖比剛才那個圖更好看，信息量也更多：它用顏色幫我們區(qū)分了組別，并添加了橫縱坐標的標簽。
讓我們添加曲線的置信區(qū)間，并增加long-rank檢驗的結果p值以及風險表格：
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曲線用來對兩個分類變量差異的可視化非常合適，但分類要是多，那就糟透了：
ggsurvplot(survfit(Surv(time, status)~nodes, data=survival::colon))
而且生存曲線另外不能可視化的是連續(xù)型變量的風險。
Cox PH回歸模型正好是處理這類問題的一把好手，它同樣內置于survival包中，語法與lm()和glm()一致。
讓我們再來用肺癌數(shù)據(jù)集看看不同性別的風險，這次使用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
結果中的exp(coef)列包含eβ1eβ1。它就是風險比率——該變量對風險率的乘數(shù)效應（對于該變量每個單位增加的）。因此，對于像性別這樣的分類變量，從男性（基線）到女性的結果大約減少約40％的危險。你也可以翻轉coef列上的符號，并采用exp(0.531)，你可以將其解釋為男性導致危險增加1.7倍，或者單位時間男性的死亡率約為女性1.7倍（女性死亡率為男性的0.588倍）。
要記住：
- HR=1: 沒有效應
- HR>1: 風險增加
- HR<1: 風險減少 （保護變量）
你還會注意到，“性別”有一個對應的p值，整個模型中也有一個p值。0.00111的p值非常接近我們在Kaplan-Meier圖上看到的p=0.00131的p值。這是因為KM曲線顯示的是對數(shù)秩檢驗的p值。你可以通過調用summary(fit)來獲得Cox模型結果。你也可以使用survdiff()直接計算log-rank測試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)建另一個模型，分析數(shù)據(jù)集中的所有變量！這向我們展示了當所有變量一起考慮時，如何影響生存。一些是非常強大的預測指標（性別，ECOG評分）。有趣的是，醫(yī)師對Karnofsky表現(xiàn)評分的評分稍高，但患者評分相同。
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ù)變量時似乎有一點重要。
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圖。但是，正如我們之前所看到的，我們不能這樣做，因為我們會為每個獨特的年齡值獲得單獨的曲線！
ggsurvplot(survfit(Surv(time, status)~age, data=lung))
你可能在這里看到的一件事是試圖將一個連續(xù)變量分成不同的組 - 三分位數(shù)，上四分位數(shù)與下四分位數(shù)，中位數(shù)分數(shù)等 - 這樣你就可以生成生存曲線圖。但是，你如何進行分組是有意義的！檢查cut的幫助。cut()接受一個連續(xù)變量和一些斷點，并從中創(chuàng)建一個分類變量。 我們來得到數(shù)據(jù)集的平均年齡，并繪制一個顯示年齡分布的直方圖。
mean(lung$age)
hist(lung$age)
ggplot(lung, aes(age)) + geom_histogram(bins=20)
現(xiàn)在，讓我們嘗試通過lung$age創(chuàng)建一個分類變量，其中0,62（平均值）和正無窮大。我們可以在這里繼續(xù)添加labels =選項來標記我們創(chuàng)建的分組，例如，“年輕”和“老”。最后，我們可以將這個結果分配給肺數(shù)據(jù)集中的一個新對象。
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)在，當我們用這個新的分類生成KM圖時會發(fā)生什么？ 看起來“老”和“年輕”患者之間的曲線存在一些差異，老年患者的生存幾率稍差。但是p=0.39時，62歲以下和62歲以上者的生存率差異不顯著。
ggsurvplot(survfit(Surv(time, status)~agecat, data=lung), pval=TRUE)
但是，如果我們選擇一個不同的切點，例如70歲，這大概是年齡分布上限的四分之一（見“分位數(shù)”）。結果現(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)
請記住，Cox回歸分析整個分布范圍內的連續(xù)變量，其中Kaplan-Meier圖上的對數(shù)秩檢驗可能會根據(jù)您對連續(xù)變量進行分類而發(fā)生變化。他們以一種不同的方式回答類似的問題：回歸模型提出的問題是“年齡對生存有什么影響？”，而對數(shù)秩檢驗和KM圖則問：“那些不到70歲和70歲以上的人有差異嗎？“。
在新的survminer 0.2.4版本中，新增了可以一次確定一個或多個連續(xù)變量最佳分割點的函數(shù)surv_cutpoint()與surv_categorize()。參閱這篇博文學習詳細的函數(shù)用法。
英文來源：http:///r-survival.html。上面講述的內容涵蓋了R生存分析的基本各個方面，具體的使用、圖形優(yōu)化還需要讀者自己琢磨。英文來源中有不少練習，有興趣的不妨試試。后續(xù)還有一些TCGA的分析實例，有空我再寫寫。
- 在醫(yī)學界，我們通常會從字面上思考生存分析 - 追蹤直至死亡的時間。 但是，它比這更普遍-生存分析模擬事件發(fā)生之前的時間（任何事件）都可以。這可以是生物有機體的死亡。但也可能是機械系統(tǒng)發(fā)生硬件故障，恢復時間，失去工作后有人失業(yè)的時間，直到成熟的番茄被放牧的鹿吃掉的時間，直到有人在車間里睡著的時間， 生存分析還包括工程可靠性理論，經濟學持續(xù)時間分析和社會學事件歷史分析。?
- 這描述了最常見的截尾類型 - 右截尾。當“開始”未知時，例如當初始診斷或暴露時間未知時，通常不會發(fā)生左側截尾。?
- 而且，按照上述定義，假定兩組之間的累積危險比率隨著時間的推移保持不變。?
- 對于這些差異，有一個類似卡方的統(tǒng)計檢驗，稱為對數(shù)秩檢驗，比較生存函數(shù)分類組。?
- Cox回歸和來自
survdiff的logrank檢驗將在大多數(shù)時間給你類似的結果。對數(shù)秩檢驗是詢問兩組患者生存曲線是否顯著不同。Cox回歸是詢問許多分類或連續(xù)變量中哪一個顯著影響生存。?
-
Surv()也可以輸入開始與截止時間，參見?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.?
|
