生存分析與R

2018年05月19日 19:55:06 走在碼農路上的醫學狗 閱讀數：4399更多

個人分類： R語言

 

版權聲明：本文為博主原創文章，未經博主允許不得轉載。 https://blog.csdn.net/weixin_40466280/article/details/80377317

生存分析與R

生存分析是將事件的結果和出現這一結果所經歷的時間結合起來分析的一類統計分析方法。不僅考慮事件是否出現，而且還考慮事件出現的時間長短，因此這類方法也被稱為事件時間分析（time-to-event analysis）。生存分析是醫學領域中一個重要的內容，在腫瘤等疾病的研究中運用十分廣泛。

1.生存分析中的重要概念

生存分析的數據資料與其它一般的數據資料有一些不同的特征： 
1. 其同時考慮生存時間和生存結局 
2. 通常存在刪失（censored）數據 
3. 生存時間通常不服從生態分布。

1.1 生存時間

生存時間（survival time）指的是從開始事件到終點事件所經歷的事件跨度。例如，急性白血病患者從發病到死亡所經歷的事件跨度，冠心病患者兩次發作之間的時間間隔等。 
注意：在進行實驗設計時，需要對起始事件、終點事件、時間單位進行明確的定義。

1.2 刪失

生存結局（status）一般分為「死亡」和刪失兩類。「死亡」指的是我們感興趣的終點事件（如白血病患者死亡、冠心病患者第二次發病）。除此之外的結局或生存結局則歸類為刪失（censoring），也稱為截尾或終檢。 
刪失的一般原因有： 
1. 研究截至日期時，感興趣終點事件仍未出現 
2. 失訪，不知道感興趣終點事件何時發生或是否會發生 
3. 因各種原因中途退出 
4. 死於其它「事件」，如交通意外或其他疾病

2 生存分析的統計學方法與R的實現

生存分析擁有着與其它分析不同的統計學方法。 
1. 描述統計：常采用Kaplan-Meier法進行分析，並繪制生存曲線；對於頻數表資料，則可以采用壽命表進行分析（屬於非參數統計方法） 
2. 比較分析：我們經常需要對不同組別的生存率進行比較分析，比如比較使用或不用某種葯物的HIV陽性患者的生存率是否不同。經常采用的log-rank檢驗以及Breslow檢驗。檢驗的零假設為：兩組或多組總體生存時間分布相同。 
3. 影響因素分析：我們可以建立生存模型來探討哪些因素影響生存時間。常用的方法有兩類，一類為半參數法：Cox比例風險模型；還有一類為參數法，主要有logistic分布法、Gompertz分布法等回歸模型。

2.1 用R繪制生存曲線

在R中進行生存分析常用的包有survival包以及survminer包。 
- survival 包提供了生存函數的建立，Cox模型的建立，以及比較分析。這個包也提供了基於基礎繪圖系統的生存曲線繪制。 
- * survminer包*提供了基於ggplot2系統的可視化，具有更加美觀的圖形，以及定制方式。

2.1.1 數據集（data set）

我們在這里使用的數據集是survival包中含有的肺癌數據集：lung。前6條數據如下：

> library(survival)
    ## 前6條數據 > head(lung) inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss 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

解釋： 
- 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 (bad=0-good=100) rated by physician Karnofsky 
- pat.karno: performance score as rated by patient 
- meal.cal: Calories consumed at meals 
- wt.loss: Weight loss in last six months

2.1.2 生存曲線的擬合

survival包中的Sruv 函數可以創建一個生存對象。

    >fit.surv <-Surv(lung$time,lung$status) > head(fit.surv) [1] 306 455 1010+ 210 883 1022+

該函數根據是否為刪失，將時間進行分類。右側帶有「+」號，代表為刪失數據。 
這時，我們可以使用survival包中的survfit函數用Kaplan-Meier法進行生存曲線的擬合。

> km<-survfit(fit.surv~1,data = lung)

同時，我們也可以將其根據年齡（age）分為兩組進行擬合：

>km_2<- survfit(fit.surv~sex,data=lung)

2.1.3 生存曲線可視化的方法：

可視化方法有兩種，一種是基於基礎繪圖系統的plot，還有一種是基於ggplot2繪圖系統的ggsurvplot。 
我們可以使用基於基礎繪圖系統的plot函數將其可視化：

    plot (km)

 
分組的生存曲線：

plot (km_2)

 
這種基於基礎繪圖系統的可視化方法較為簡陋，可以修改的也參數也較少。目前在R語言中，可視化功能極其強大的是ggplot2系統。survminer包提供了基於ggplot2系統的可視化函數：ggsurvplot

    library(survminer) ggsurvplot (km)

ggsurvplot (km_2)

 
對比可以看出，survminer包繪制的生存曲線更加美觀。同樣，我們還可以對圖形進行顏色的改變，圖形大小的改變，增加圖標以及圖例位置的更換等：

    ggsurvplot(km_2,
               legend = "bottom", #將圖例移動到下方 legend.title = "Sex",#改變圖例名稱 legend.labs = c("Male", "Female"), linetype = "strata"# 改變線條類型 )

修改后的圖形： 

3.比較分析

一般情況下，我們都需要比較分析兩組的生存時間分布是否不同。在survival包中，有一個survdiff的函數可以進行long-rank檢驗

    > survdiff(fit.surv~sex,data = lung,
    rho = 0 # rho = 0 表示使用long-rank檢驗或者Mantel-Haenszel 檢驗) Call: survdiff(formula = fit.surv ~ sex, data = lung, rho = 0) 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 > 

此外，可以使用survminer包中的ggsurvplot函數中的pval=TRUE參數，在生存曲線中添加P值：

    ggsurvplot(km_2, main = "Survival curve", pval=TRUE #添加P值 )

4.Cox 回歸模型

4.1 Cox 模型的建立

我們還是利用第二部分所用的肺癌數據集（lung）進行Cox回歸模型的建立，這一次，我們感興趣的點主要是年齡、體重減輕以及性別是否會影響肺癌的生存時間：

    > library(survival)
    > res.cox<-coxph(Surv(time,status)~age+ph.ecog+wt.loss,data = lung) > res.cox Call: coxph(formula = Surv(time, status) ~ age + ph.ecog + wt.loss, data = lung) coef exp(coef) se(coef) z p age 0.01347 1.01356 0.00974 1.38 0.16659 ph.ecog 0.47222 1.60356 0.12771 3.70 0.00022 wt.loss -0.00717 0.99285 0.00663 -1.08 0.27921 Likelihood ratio test=19 on 3 df, p=0.000269 n= 213, number of events= 151 (15 observations deleted due to missingness) 

在這個模型中我們可以看到ECOG評分的P值\<0.01，認為在這個模型中能夠顯著影響肺癌的生存分析，OR=1.6。

4.2 模型診斷——PH檢驗

在構建Cox模型之后，我們需要對這個模型進行PH檢驗。如果無法通過PH檢驗，我們需要對上述的Cox模型進行修改。 
在survival包中，函數cox.zph可進行PH檢驗：

cox.zph(res.cox) rho chisq p age -0.03663 0.22364 0.636 ph.ecog -0.08018 1.20971 0.271 wt.loss -0.00191 0.00064 0.980 GLOBAL NA 1.95904 0.581 > 

可以看到P 值都>0.01，說明該模型能夠通過PH檢驗。 
使用survminer包中的ggcoxzph函數還可以將其進行可視化：

    > temp<-cox.zph(res.cox) > ggcoxzph(temp)

如果無法通過PH檢驗，可以進行分層或使用其他的檢驗方法，如參數檢驗等。