基本概念
生存分析：從字面上就是讓我們分析事件發(fā)生的速率，研究各個因素與生存時間有無關(guān)系及關(guān)聯(lián)程度大小。主要描述3個內(nèi)容：
1.研究生存過程
研究生存時間的分布特點，估計生存率及標準誤、繪制生存曲線、繪制生存曲線。常用K-M法（乘積極限法）和壽命法。生存分析為單因素分析（兩個或多個水平），用中位生存時間即50%生存率對應的生存時間表示生存時間的平均水平。
2.比較生存過程
獲得生存率及其標準誤的估計值后，可進行兩組或多組生存曲線的比較。常用方法有對數(shù)秩檢驗（log-rank檢驗），若曲線交叉，可能存在混雜因素。
3.影響生存時間的因素分析
常用的多因素生存分析方法：cox比例風險回歸模型。
起始事件
終點事件
觀察時間
生存時間：用t表示
完全數(shù)據(jù)
截尾數(shù)據(jù)（刪失值）：觀察時間不是由于終點事件結(jié)束的，可能是失訪、死于非研究因素、觀察結(jié)束而對象仍存活的結(jié)局。
風險是特定時間點t的瞬時事件發(fā)生（死亡）率。生存分析并不認為隨著時間的推移危害是不變的。累積風險是直到時間t為止經(jīng)歷的總風險。生存函數(shù)是個體在時間t之前存在的概率（或者，不發(fā)生感興趣事件的概率）。這是事件（例如，死亡）尚未發(fā)生的可能性?？雌饋硐襁@樣，其中T是死亡時間，Pr(T>t)是死亡時間大于某個時間t的概率。S是概率，所以0≤S(t)≤1，因為生存期總是正值（T≥0）。
S(t)=Pr(T>t)。
Kaplan-Meier曲線描述了生存函數(shù)。呈現(xiàn)了不同組間患者生存概率隨時間的變化，能很好地描述生存過程。這是一個階梯函數(shù)，說明隨著時間的推移累計生存概率。曲線在沒有事件發(fā)生的時間段內(nèi)是水平的，然后垂直下降，對應于每次發(fā)生事件時生存函數(shù)的變化。截尾是一種生存分析特有的缺失數(shù)據(jù)問題。 當我們在研究結(jié)束時跟蹤樣本/主題并且事件從未發(fā)生時會發(fā)生這種情況。這也可能是由于樣本/受試者因死亡以外的原因而退出研究或其他一些失訪導致的。樣本數(shù)據(jù)發(fā)生截尾是因為你只知道這個人存活到失去跟蹤為止，但你不知道任何關(guān)于之后他的生存狀態(tài)。
比例風險假設(shè)：生存分析的主要目標是比較不同組別（例如白血病患者與正常對照組）的生存功能。如果兩組人都死亡，那么兩條生存曲線都將結(jié)束于0％，但是其中一組的平均存活時間可能比另一組長。生存分析通過比較觀察期間不同時間的風險來做到這一點。生存分析并不假設(shè)危害是恒定的，但假定組間危害的比率隨著時間的推移是恒定的。
比例風險回歸也稱為Cox回歸，是評估不同變量對生存率影響的最常用方法。
Cox PH模型
Kaplan-Meier曲線適用于觀察兩個分類組之間生存率的差異，但對于評估諸如年齡，基因表達，白細胞計數(shù)等定量變量的影響，它們不起作用。Cox PH回歸可評估分類變量和連續(xù)變量的影響，并且可以一次模擬多個變量的影響。
Cox PH回歸模型將時間t處的風險自然對數(shù)表示為h(t)，作為基線危險（h0(t)）的函數(shù)（所有暴露變量為0的個體的風險）和多個暴露變量x1x1，x1x1，......，xpxp。 Cox PH模型的形式是：
log(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內(nèi)保持不變（因此名稱比例風險回歸）。ββ值是根據(jù)模型估計的回歸系數(shù)，并表示相應預測變量中每單位增加的log(Hazard,Ratio)log(Hazard,Ratio)。危害比的解釋取決于預測變量的測量尺度，但簡單地說，正系數(shù)表示較差的存活率，負系數(shù)表示所討論的變量存活率較高。
所以一般做生存分析，可以用KM（Kaplan-Meier）方法估計生存率，做生存曲線，然后可以根據(jù)分組檢驗一下多組間生存曲線是否有顯著的差異，最后用Cox風險比例模型來研究下某個因素對生存的影響
R做生存分析：
用survival包做生存分析、survminer的總結(jié)和可視化生存分析結(jié)果
Log-Rank檢驗比較生存曲線：survdiff（），Log-Rank test是無參數(shù)檢驗，近似于卡方檢驗，零假設(shè)是組間沒有差異
對數(shù)秩檢驗是比較兩條或更多條生存曲線的最廣泛使用的方法。零假設(shè)是兩組在生存期間沒有差異。對數(shù)秩檢驗是一個非參數(shù)檢驗，不存在關(guān)于生存分布的假設(shè)。從本質(zhì)上講，對數(shù)秩檢驗將觀察到的每組事件數(shù)量與假設(shè)假設(shè)為真（即生存曲線是否相同）所預期的數(shù)量進行比較。對數(shù)秩的統(tǒng)計近似地分布為卡方檢驗統(tǒng)計量。
存活包中的函數(shù)survdiff（）可以用來計算比較兩個或更多生存曲線的log-rank測試。
surv_diff<-survdiff(Surv(time,status)~sex,data=lung)
surv_diff
#Call:survdiff(formula = Surv(time, status) ~ sex, data = lung)N Observed Expected (O-E)^2/E (O-E)^2/Vsex=1 138 112 91.6 4.55 10.3sex=2 90 53 73.4 5.68 10.3Chisq= 10.3 on 1 degrees of freedom, p= 0.0013
若可視化圖可看出性別組間生存曲線是有差異的，也從Log-Rank test的P值0.0013也可得出性別組有顯著差異.
可視化生存曲線
我們將使用函數(shù)ggsurvplot（）（在SurvminerR軟件包中）來生成兩組受試者的生存曲線。
也可以顯示：
使用參數(shù)conf.int = TRUE的幸存函數(shù)的95％置信限。
使用期權(quán)risk.table的風險個體的數(shù)量和/或百分比。risk.table允許的值包括：
指定是否顯示風險表的TRUE或FALSE。默認值是FALSE。
“絕對”或“百分比”：分別顯示風險對象的絕對數(shù)量和百分比。使用“abs_pct”來顯示絕對數(shù)字和百分比。
Log-Rank測試的p值使用pval = TRUE比較組。
在使用參數(shù)surv.median.line的中值生存中的水平/垂直線。
好啦，接下來做用R做一遍生存分析
核心的分析函數(shù)都在survival包里，我們這里使用dplyr包，然后用survminer包進行繪圖。
library(survival)
library(dplyr)
library(survminer)
核心函數(shù)包括：
-Surv()：創(chuàng)建一個生存對象
-survfit()：使用公式或已構(gòu)建的Cox模型擬合生存曲線
-coxph()：擬合Cox比例風險回歸模型
其他我們可能會用到的函數(shù)：
-
cox.zph()：檢驗一個Cox回歸模型的比例風險假設(shè) -
survdiff()：用log-rank/Mantel-Haenszel檢驗檢驗生存差異5 -
Surv()創(chuàng)建響應變量，典型用法是使用事件時間，6以及事件是否發(fā)生（即死亡與截尾）。 -
survfit()創(chuàng)建一條生存曲線，然后可以顯示或繪圖。 -
coxph()實現(xiàn)回歸分析，并且模型以與常規(guī)線性模型中相同的方式指定，但使用coxph()函數(shù)。
我們將使用內(nèi)置的肺癌數(shù)據(jù)集7學習使用survival包。
library(survival)
?lung#獲取數(shù)據(jù)集信息，可在右手邊查看
Format
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 as rated by the physician. 0=asymptomatic, 1= symptomatic but completely ambulatory, 2= in bed <50% of the day, 3= in bed > 50% of the day but not bedbound, 4 = bedbound
ph.karno: Karnofsky performance score (bad=0-good=100) 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)#我覺得data(lung)也可以
lung
## A tibble: 228 x 10
inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 3 306 2 74 1 1 90 100 1175 NA
2 3 455 2 68 1 0 90 90 1225 15
3 3 1010 1 56 1 0 90 90 NA 15
4 5 210 2 57 1 1 90 60 1150 11
5 1 883 2 60 1 0 100 90 NA 0
6 12 1022 1 74 1 1 50 80 513 0
7 7 310 2 68 2 2 70 60 384 10
8 11 361 2 71 2 2 60 80 538 1
9 1 218 2 53 1 1 70 80 825 16
10 7 166 2 61 1 2 70 70 271 34
# ... with 218 more rows
生存曲線
1.構(gòu)建生存對象：
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 654 728 71 567 144 613 707 61 88 301 81 624 371 394 520
[27] 574 118 390 12 473 26 533 107 53 122 814 965+ 93 731 460 153 433 145 583 95 303 519 643 765 735 189
[53] 53 246 689 65 5 132 687 345 444 223 175 60 163 65 208 821+ 428 230 840+ 305 11 132 226 426 705 363
[79] 11 176 791 95 196+ 167 806+ 284 641 147 740+ 163 655 239 88 245 588+ 30 179 310 477 166 559+ 450 364 107
[105] 177 156 529+ 11 429 351 15 181 283 201 524 13 212 524 288 363 442 199 550 54 558 207 92 60 551+ 543+
[131] 293 202 353 511+ 267 511+ 371 387 457 337 201 404+ 222 62 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 181 285 301+ 348 197 382+ 303+ 296+ 180 186 145 269+ 300+ 284+
[183] 350 272+ 292+ 332+ 285 259+ 110 286 270 81 131 225+ 269 225+ 243+ 279+ 276+ 135 79 59 240+ 202+ 235+ 105 224+ 239
[209] 237+ 173+ 252+ 221+ 185+ 92+ 13 222+ 192+ 183 211+ 175+ 197+ 203+ 116 188+ 191+ 105+ 174+ 177+
2.用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
# 可以把前面操作即構(gòu)建對象、擬合曲線一步完成
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ù)查看模型匯總結(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
## [到達getOption("max.print") -- 略過111行]]
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
summary()函數(shù)中可以設(shè)定時間參數(shù)用來選定一個時間區(qū)間，我們可以以此比對男生是不是比女生有更高的風險：
summary(sfit, times=seq(0, 1000, 100))#取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)
survminer的包提供的一個叫g(shù)gsurvplot()的函數(shù)可以幫助我們更簡單地做出可以發(fā)表的生存曲線,
library(survminer)
ggsurvplot(sfit)
這個圖比剛才那個圖更好看，信息量也更多：它用顏色幫我們區(qū)分了組別，并添加了橫縱坐標的標簽。
讓我們添加曲線的置信區(qū)間，并增加long-rank檢驗的結(jié)果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=.3)
Cox回歸模型
Cox PH回歸模型可以可視化連續(xù)型變量的風險，它同樣內(nèi)置于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.5310 0.5880 0.1672 -3.176 0.00149
#Likelihood ratio test=10.63 on 1 df, p=0.001111
#n= 228, number of events= 165
結(jié)果中的exp(coef)列包含eβ1eβ1。它就是風險比率——該變量對風險率的乘數(shù)效應（對于該變量每個單位增加的）。因此，對于像性別這樣的分類變量，從男性（基線）到女性的結(jié)果大約減少約40％的危險。你也可以翻轉(zhuǎn)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值。你可以通過調(diào)用summary(fit)來獲得Cox模型結(jié)果。你也可以使用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.5310 0.5880 0.1672 -3.176 0.00149 **
#---
#Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# exp(coef) exp(-coef) lower .95 upper .95
#sex 0.588 1.701 0.4237 0.816
#Concordance= 0.579 (se = 0.021 )
#Likelihood ratio test= 10.63 on 1 df, p=0.001
#Wald test = 10.09 on 1 df, p=0.001
#Score (logrank) test = 10.33 on 1 df, p=0.001
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.001
本來p應該是0.0013的，但是不知道是錯在哪了，還是因為取約數(shù)問題？
創(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
#all:
#oxph(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
#ex -5.509e-01 5.765e-01 2.008e-01 -2.743 0.00609
#ge 1.065e-02 1.011e+00 1.161e-02 0.917 0.35906
#h.ecog 7.342e-01 2.084e+00 2.233e-01 3.288 0.00101
#h.karno 2.246e-02 1.023e+00 1.124e-02 1.998 0.04574
#at.karno -1.242e-02 9.877e-01 8.054e-03 -1.542 0.12316
#eal.cal 3.329e-05 1.000e+00 2.595e-04 0.128 0.89791
#t.loss -1.433e-02 9.858e-01 7.771e-03 -1.844 0.06518
#ikelihood ratio test=28.33 on 7 df, p=0.0001918
#= 168, number of events= 121
# (60 observations deleted due to missingness)
分類畫KM曲線
繼續(xù)看年齡的Cox模型。看起來年齡在模擬為連續(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.018720 1.018897 0.009199 2.035 0.0419
#Likelihood ratio test=4.24 on 1 df, p=0.03946
#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)#求年齡的平均數(shù)
hist(lung$age)#畫直方圖
ggplot(lung, aes(age)) + geom_histogram(bins=20)
現(xiàn)在，嘗試通過lung$age創(chuàng)建一個分類變量，其中0,62（平均值）和正無窮大。我們可以在這里繼續(xù)添加labels =選項來標記我們創(chuàng)建的分組，例如，“年輕”和“老”。最后，我們可以將這個結(jié)果分配給肺數(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] (62,Inf] (0,62] (0,62] (0,62] (62,Inf] (62,Inf] (0,62] (0,62] (62,Inf] (62,Inf]
# [18] (62,Inf] (0,62] (0,62] (62,Inf] (0,62] (0,62] (0,62] (62,Inf] (62,Inf] (0,62] (62,Inf] (0,62] (62,Inf] (62,Inf] (62,Inf] (0,62] (0,62]
# [35] (0,62] (0,62] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (0,62] (0,62] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (0,62] (62,Inf]
# [52] (62,Inf] (62,Inf] (0,62] (0,62] (0,62] (62,Inf] (0,62] (0,62] (62,Inf] (62,Inf] (0,62] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (62,Inf]
# [69] (62,Inf] (62,Inf] (62,Inf] (0,62] (62,Inf] (0,62] (0,62] (62,Inf] (0,62] (0,62] (62,Inf] (62,Inf] (0,62] (0,62] (0,62] (0,62] (0,62]
#[86] (62,Inf] (0,62] (0,62] (0,62] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (0,62] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (0,62] (62,Inf]
#[103] (0,62] (62,Inf] (0,62] (62,Inf] (0,62] (62,Inf] (0,62] (62,Inf] (62,Inf] (0,62] (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] (0,62] (62,Inf] (62,Inf] (0,62] (0,62] (0,62] (0,62] (0,62] (62,Inf] (62,Inf]
#[137] (0,62] (0,62] (0,62] (0,62] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (0,62] (0,62] (62,Inf] (0,62] (62,Inf] (0,62] (62,Inf] (0,62] (0,62]
#[154] (0,62] (0,62] (62,Inf] (62,Inf] (0,62] (62,Inf] (0,62] (0,62] (0,62] (62,Inf] (62,Inf] (62,Inf] (0,62] (0,62] (0,62] (0,62] (62,Inf]
#[171] (0,62] (0,62] (0,62] (0,62] (0,62] (0,62] (0,62] (0,62] (0,62] (62,Inf] (0,62] (0,62] (62,Inf] (62,Inf] (0,62] (0,62] (62,Inf]
#[188] (0,62] (62,Inf] (0,62] (62,Inf] (0,62] (0,62] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (0,62] (0,62] (62,Inf] (62,Inf] (62,Inf] (0,62]
#[205] (62,Inf] (0,62] (0,62] (0,62] (62,Inf] (0,62] (0,62] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (62,Inf] (0,62] (62,Inf] (62,Inf] (0,62] (62,Inf]
#[222] (62,Inf] (62,Inf] (62,Inf] (0,62] (62,Inf] (62,Inf] (0,62]
#Levels: (0,62] (62,Inf]
cut(lung$age, breaks=c(0, 62, Inf), labels=c("young", "old"))#區(qū)分年老
# [1] old old young young young old old old young young young old old young young old old old young young old young young young old old
# [27] young old young old old old young young young young old old old old old old young young old old old old old young old old
# [53] old young young young old young young old old young old old old old old old old old old young old young young old young young
# [79] old old young young young young young old young young young old old old old young old old old old old old young old young old
#[105] young old young old young old old young old old young old young old old old old young old old old old young old old young
#[131] young young young young old old young young young young old old old old young young old young old young old young young young young old
#[157] old young old young young young old old old young young young young old young young young young young young young young young old young young
#[183] old old young young old young old young old young young old old old old old young young old old old young old young young young
#[209] old young young old old old old old young old old young old old old old young old old young
#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 wt.loss agecat
# <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <fct>
#1 3 306 2 74 1 1 90 100 1175 NA old
#2 3 455 2 68 1 0 90 90 1225 15 old
#3 3 1010 1 56 1 0 90 90 NA 15 young
#4 5 210 2 57 1 1 90 60 1150 11 young
#5 1 883 2 60 1 0 100 90 NA 0 young
#6 12 1022 1 74 1 1 50 80 513 0 old
現(xiàn)在，當我們用這個新的分類生成KM圖時會發(fā)生什么？ 看起來“老”和“年輕”患者之間的曲線存在一些差異，老年患者的生存幾率稍差。但是p=0.39時，62歲以下和62歲以上者的生存率差異不顯著。
ggsurvplot(survfit(Surv(time, status)~agecat, data=lung), pval=TRUE)
但是，如果我們選擇一個不同的切點，例如70歲，這大概是年齡分布上限的四分之一（見“分位數(shù)”）。結(jié)果變得有意義了，兩組有顯著的差異性。
# 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回歸分析整個分布范圍內(nèi)的連續(xù)變量，其中Kaplan-Meier圖上的對數(shù)秩檢驗可能會根據(jù)您對連續(xù)變量進行分類而發(fā)生變化。他們以一種不同的方式回答類似的問題：回歸模型提出的問題是“年齡對生存有什么影響？”，而對數(shù)秩檢驗和KM圖則問：“那些不到70歲和70歲以上的人有差異嗎？”。
參考文章：
1.[王詩翔 【r<-統(tǒng)計|繪圖】使用R進行生存分析——一文打盡]（http://www.itdecent.cn/p/4ad9ba730719?utm_campaign=haruki&utm_content=note&utm_medium=reader_share&utm_source=qq）
https://www.mediecogroup.com/method_topic_article_detail/173/
https://zhuanlan.zhihu.com/p/34802509
2.R語言生存分析入門
https://www.bioinfo-scrounger.com/archives/647/
