基本概念
生存分析：從字面上就是讓我們分析事件發(fā)生的速率器腋，研究各個因素與生存時間有無關(guān)系及關(guān)聯(lián)程度大小脆诉。主要描述3個內(nèi)容：
1.研究生存過程
研究生存時間的分布特點雾棺，估計生存率及標(biāo)準(zhǔn)誤卸夕、繪制生存曲線柱告、繪制生存曲線截驮。常用K-M法（乘積極限法）和壽命法。生存分析為單因素分析（兩個或多個水平）际度，用中位生存時間即50%生存率對應(yīng)的生存時間表示生存時間的平均水平葵袭。
2.比較生存過程
獲得生存率及其標(biāo)準(zhǔn)誤的估計值后，可進行兩組或多組生存曲線的比較乖菱。常用方法有對數(shù)秩檢驗（log-rank檢驗）坡锡，若曲線交叉，可能存在混雜因素窒所。
3.影響生存時間的因素分析
常用的多因素生存分析方法：cox比例風(fēng)險回歸模型鹉勒。
起始事件
終點事件
觀察時間
生存時間：用t表示
完全數(shù)據(jù)
截尾數(shù)據(jù)（刪失值）：觀察時間不是由于終點事件結(jié)束的，可能是失訪吵取、死于非研究因素禽额、觀察結(jié)束而對象仍存活的結(jié)局。
風(fēng)險是特定時間點t的瞬時事件發(fā)生（死亡）率皮官。生存分析并不認為隨著時間的推移危害是不變的脯倒。累積風(fēng)險是直到時間t為止經(jīng)歷的總風(fē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)是水平的逼龟，然后垂直下降，對應(yīng)于每次發(fā)生事件時生存函數(shù)的變化追葡。截尾是一種生存分析特有的缺失數(shù)據(jù)問題腺律。 當(dāng)我們在研究結(jié)束時跟蹤樣本/主題并且事件從未發(fā)生時會發(fā)生這種情況。這也可能是由于樣本/受試者因死亡以外的原因而退出研究或其他一些失訪導(dǎo)致的宜肉。樣本數(shù)據(jù)發(fā)生截尾是因為你只知道這個人存活到失去跟蹤為止匀钧，但你不知道任何關(guān)于之后他的生存狀態(tài)。
比例風(fēng)險假設(shè)：生存分析的主要目標(biāo)是比較不同組別（例如白血病患者與正常對照組）的生存功能谬返。如果兩組人都死亡之斯，那么兩條生存曲線都將結(jié)束于0％，但是其中一組的平均存活時間可能比另一組長遣铝。生存分析通過比較觀察期間不同時間的風(fēng)險來做到這一點佑刷。生存分析并不假設(shè)危害是恒定的，但假定組間危害的比率隨著時間的推移是恒定的酿炸。
比例風(fēng)險回歸也稱為Cox回歸瘫絮，是評估不同變量對生存率影響的最常用方法。
Cox PH模型
Kaplan-Meier曲線適用于觀察兩個分類組之間生存率的差異梁沧，但對于評估諸如年齡檀何，基因表達，白細胞計數(shù)等定量變量的影響廷支，它們不起作用频鉴。Cox PH回歸可評估分類變量和連續(xù)變量的影響，并且可以一次模擬多個變量的影響恋拍。
Cox PH回歸模型將時間t處的風(fēng)險自然對數(shù)表示為h(t)垛孔，作為基線危險（h0(t)）的函數(shù)（所有暴露變量為0的個體的風(fēng)險）和多個暴露變量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
重新排列該等式可以估計風(fēng)險比率默怨，比較時間t暴露對于未暴露個體：
HR(t)=h1(t)h0(t)=eβ1HR(t)=h1(t)h0(t)=eβ1
該模型顯示風(fēng)險比為eβ1eβ1讯榕，并且在時間t內(nèi)保持不變（因此名稱比例風(fēng)險回歸）。ββ值是根據(jù)模型估計的回歸系數(shù)，并表示相應(yīng)預(yù)測變量中每單位增加的log(Hazard,Ratio)log(Hazard,Ratio)愚屁。危害比的解釋取決于預(yù)測變量的測量尺度济竹，但簡單地說，正系數(shù)表示較差的存活率霎槐，負系數(shù)表示所討論的變量存活率較高送浊。
所以一般做生存分析，可以用KM（Kaplan-Meier）方法估計生存率丘跌，做生存曲線袭景，然后可以根據(jù)分組檢驗一下多組間生存曲線是否有顯著的差異，最后用Cox風(fēng)險比例模型來研究下某個因素對生存的影響
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è)為真（即生存曲線是否相同）所預(yù)期的數(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的風(fēng)險個體的數(shù)量和/或百分比一睁。risk.table允許的值包括：
指定是否顯示風(fēng)險表的TRUE或FALSE。默認值是FALSE佃却。
“絕對”或“百分比”：分別顯示風(fēng)險對象的絕對數(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比例風(fēng)險回歸模型
其他我們可能會用到的函數(shù)：
-
cox.zph()：檢驗一個Cox回歸模型的比例風(fēng)險假設(shè) -
survdiff()：用log-rank/Mantel-Haenszel檢驗檢驗生存差異5 -
Surv()創(chuàng)建響應(yīng)變量赦邻，典型用法是使用事件時間髓棋，6以及事件是否發(fā)生（即死亡與截尾）。 -
survfit()創(chuàng)建一條生存曲線深纲，然后可以顯示或繪圖仲锄。 -
coxph()實現(xiàn)回歸分析劲妙，并且模型以與常規(guī)線性模型中相同的方式指定，但使用coxph()函數(shù)儒喊。
我們將使用內(nèi)置的肺癌數(shù)據(jù)集7學(xué)習(xí)使用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ū)間，我們可以以此比對男生是不是比女生有更高的風(fēng)險：
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ū)分了組別，并添加了橫縱坐標(biāo)的標(biāo)簽绑警。
讓我們添加曲線的置信區(qū)間求泰，并增加long-rank檢驗的結(jié)果p值以及風(fēng)險表格：
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ù)型變量的風(fēng)險，它同樣內(nèi)置于survival包中计盒，語法與lm()和glm()一致渴频。
讓我們再來用肺癌數(shù)據(jù)集看看不同性別的風(fēng)險，這次使用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卜朗。它就是風(fēng)險比率——該變量對風(fēng)險率的乘數(shù)效應(yīng)（對于該變量每個單位增加的）。因此咕村，對于像性別這樣的分類變量场钉，從男性（基線）到女性的結(jié)果大約減少約40％的危險。你也可以翻轉(zhuǎn)coef列上的符號培廓，并采用exp(0.531)惹悄，你可以將其解釋為男性導(dǎo)致危險增加1.7倍，或者單位時間男性的死亡率約為女性1.7倍（女性死亡率為男性的0.588倍）肩钠。
注意
HR=1: 沒有效應(yīng)
HR>1: 風(fēng)險增加
HR<1: 風(fēng)險減少 （保護變量）
我們還注意到泣港，“性別”有一個對應(yīng)的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應(yīng)該是0.0013的，但是不知道是錯在哪了，還是因為取約數(shù)問題箫柳？
創(chuàng)建另一個模型手形，分析數(shù)據(jù)集中的所有變量！這向我們展示了當(dāng)所有變量一起考慮時悯恍，如何影響生存库糠。一些是非常強大的預(yù)測指標(biāo)（性別，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 =選項來標(biāo)記我們創(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)在，當(dāng)我們用這個新的分類生成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.reibang.com/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/
