生存分析與R

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

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

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

1.1 生存時(shí)間

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

1.2 刪失

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

2 生存分析的統(tǒng)計(jì)學(xué)方法與R的實(shí)現(xiàn)

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

2.1 用R繪制生存曲線

在R中進(jìn)行生存分析常用的包有survival包以及survminer包。?
-?survival 包提供了生存函數(shù)的建立，Cox模型的建立，以及比較分析。這個(gè)包也提供了基于基礎(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

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10

解釋：?
- 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)建一個(gè)生存對象。

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

• 1
• 2
• 3

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

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

• 1

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

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

• 1

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

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

plot (km)

• 1

?
分組的生存曲線：

plot (km_2)

• 1

?
這種基于基礎(chǔ)繪圖系統(tǒng)的可視化方法較為簡陋，可以修改的也參數(shù)也較少。目前在R語言中，可視化功能極其強(qiáng)大的是ggplot2系統(tǒng)。survminer包提供了基于ggplot2系統(tǒng)的可視化函數(shù)：ggsurvplot

library(survminer)ggsurvplot (km)

• 1
• 2

ggsurvplot (km_2)

• 1

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

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

• 1
• 2
• 3
• 4
• 5
• 6

修改后的圖形：?

3.比較分析

一般情況下，我們都需要比較分析兩組的生存時(shí)間分布是否不同。在survival包中，有一個(gè)survdiff的函數(shù)可以進(jìn)行l(wèi)ong-rank檢驗(yàn)

> survdiff(fit.surv~sex,data = lung,rho = 0 # rho = 0 表示使用long-rank檢驗(yàn)或者M(jìn)antel-Haenszel 檢驗(yàn))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 >

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11

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

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

• 1
• 2
• 3

4.Cox 回歸模型

4.1 Cox 模型的建立

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

> library(survival)> res.cox<-coxph(Surv(time,status)~age+ph.ecog+wt.loss,data = lung)> res.coxCall: 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)

• 1
• 2
• 3
• 4
• 5
• 6
• 7
• 8
• 9
• 10
• 11
• 12
• 13
• 14
• 15
• 16

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

4.2 模型診斷——PH檢驗(yàn)

在構(gòu)建Cox模型之后，我們需要對這個(gè)模型進(jìn)行PH檢驗(yàn)。如果無法通過PH檢驗(yàn)，我們需要對上述的Cox模型進(jìn)行修改。?
在survival包中，函數(shù)cox.zph可進(jìn)行PH檢驗(yàn)：

cox.zph(res.cox)rho chisq page -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 >

• 1
• 2
• 3
• 4
• 5
• 6
• 7

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

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

• 1
• 2

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

轉(zhuǎn)載于:https://www.cnblogs.com/nkwy2012/p/9241694.html

總結(jié)

以上是生活随笔為你收集整理的生存分析与R--转载的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。