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

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

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

1.1 生存時間

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

1.2 刪失

生存結(jié)局（status）一般分為「死亡」和刪失兩類?！杆劳觥怪傅氖俏覀兏信d趣的終點事件（如白血病患者死亡,、冠心病患者第二次發(fā)病）,。除此之外的結(jié)局或生存結(jié)局則歸類為刪失（censoring）,，也稱為截尾或終檢。
刪失的一般原因有：
1. 研究截至日期時,，感興趣終點事件仍未出現(xiàn)
2. 失訪,，不知道感興趣終點事件何時發(fā)生或是否會發(fā)生
3. 因各種原因中途退出
4. 死于其它「事件」，如交通意外或其他疾病

2 生存分析的統(tǒng)計學方法與R的實現(xiàn)

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

2.1 用R繪制生存曲線

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

2.1.1 數(shù)據(jù)集（data set）

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

> library(survival)
    ## 前6條數(shù)據(jù)
> 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 函數(shù)可以創(chuàng)建一個生存對象。

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

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

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

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

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

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

可視化方法有兩種,，一種是基于基礎(chǔ)繪圖系統(tǒng)的plot,，還有一種是基于ggplot2繪圖系統(tǒng)的ggsurvplot。
我們可以使用基于基礎(chǔ)繪圖系統(tǒng)的plot函數(shù)將其可視化：

    plot (km)

分組的生存曲線：

plot (km_2)

這種基于基礎(chǔ)繪圖系統(tǒng)的可視化方法較為簡陋,，可以修改的也參數(shù)也較少,。目前在R語言中，可視化功能極其強大的是ggplot2系統(tǒng),。survminer包提供了基于ggplot2系統(tǒng)的可視化函數(shù)：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的函數(shù)可以進行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函數(shù)中的pval=TRUE參數(shù),，在生存曲線中添加P值：

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

4.Cox 回歸模型

4.1 Cox 模型的建立

我們還是利用第二部分所用的肺癌數(shù)據(jù)集（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檢驗

在構(gòu)建Cox模型之后,，我們需要對這個模型進行PH檢驗,。如果無法通過PH檢驗，我們需要對上述的Cox模型進行修改,。
在survival包中,，函數(shù)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函數(shù)還可以將其進行可視化：

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

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