R語言Cox與KM生存結(jié)果分析

這篇“R語言Cox與KM生存結(jié)果分析”文章的知識(shí)點(diǎn)大部分人都不太理解，所以小編給大家總結(jié)了以下內(nèi)容，內(nèi)容詳細(xì)，步驟清晰，具有一定的借鑒價(jià)值，希望大家閱讀完這篇文章能有所收獲，下面我們一起來看看這篇“R語言Cox與KM生存結(jié)果分析”文章吧。

Cox與KM生存分析及結(jié)果解讀，查找與預(yù)后相關(guān)的biomarker

生存分析：針對(duì)于慢性病(癌癥)，因?yàn)槠錈o法在短時(shí)間內(nèi)判斷預(yù)后，不宜采用治愈率和病死率等指標(biāo)，而是需要對(duì)患者進(jìn)行隨訪，分析一定時(shí)間后患者生存或死亡的情況，這樣將事件的結(jié)果和出現(xiàn)這一結(jié)果所經(jīng)歷的時(shí)間結(jié)合起來分析的統(tǒng)計(jì)方法稱為生存分析。

生存分析KM法與Cox法異同介紹

KM 方法即Kaplan-Meier survival estimate是一種無參數(shù)方法(non-parametric)來從觀察的生存時(shí)間來估計(jì)生存概率的方法。KM生存分析模型，是單變量分析(univariable analysis)，在做單變量分析時(shí)，模型只描述了該單變量和生存之間的關(guān)系而忽略其他變量的影響。（為什么要考慮multi-variables？比如在比較兩組病人擁有和不擁有某種基因型對(duì)生存率的影響，但是其中一組的患者年齡較大，所以生存率可能受到基因型 或/和 年齡的共同影響）。同時(shí)，Kaplan-Meier方法只能針對(duì)分類變量（治療A vs 治療B，男 vs 女），不能分析連續(xù)變量對(duì)生存造成的影響。為了解決上述兩種問題，Cox比例風(fēng)險(xiǎn)回歸模型（Cox proportional hazards regression model）就被提了出來。因此，Cox比例風(fēng)險(xiǎn)回歸模型，可以分析連續(xù)變量對(duì)生存造成的影響，也可以多變量分析對(duì)生存的影響。

R語言中的生存分析R包介紹

R是數(shù)據(jù)分析常用的軟件之一，通過各種功能強(qiáng)大的R包，可以簡(jiǎn)單方便的實(shí)現(xiàn)各種分析。在R語言中，能夠進(jìn)行生存分析的R包很多，survival和survminer是其中最基本的兩個(gè)，survival負(fù)責(zé)分析，survimner負(fù)責(zé)可視化，二者相結(jié)合，可以輕松實(shí)現(xiàn)生存分析。

R語言中的KM模型分析

1. 準(zhǔn)備生存數(shù)據(jù)

對(duì)于每個(gè)個(gè)體而言，其生存數(shù)據(jù)會(huì)出現(xiàn)兩種情況，第一種是觀測(cè)到生存時(shí)間，通常用1表示，第二種則是刪失。通常用0表示。survival自帶了一個(gè)測(cè)試數(shù)據(jù)lung

每一行代表一個(gè)樣本，time表示生存時(shí)間，status表示刪失情況，這里只有1和2兩種取值，默認(rèn)排序后的第一個(gè)level對(duì)應(yīng)的值為刪失，這里則為1表示刪失。其他列為樣本對(duì)應(yīng)的性別，年齡等基本信息。

2. 進(jìn)行生存分析

這里根據(jù)性別這個(gè)二分類變量，采用KM算法來估計(jì)生存曲線，代碼如下(summary結(jié)果只顯示部分）

> library("survival")
> library("survminer")
> kmfit<-survfit(Surv(time, status) ~ sex, data =  lung)
> summary(kmfit)
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

從kmfit中summary可以看到已經(jīng)包含了每個(gè)時(shí)間點(diǎn)的生存概率，刪失等信息，通過這些信息，完全可以自己寫代碼來畫圖。為了方便，我們直接采用survminer中的函數(shù)來進(jìn)行可視化。

3. 分析結(jié)果的可視化

最基本的可視化方式如下

library("survminer")
ggsurvplot(kmfit,
              pval = TRUE, conf.int = TRUE,
              risk.table = TRUE,
              risk.table.col = "strata",
              linetype = "strata",
              surv.median.line = "hv",
              ggtheme = theme_bw(),
              palette = c("#E7B800", "#2E9FDF"))

語言中的Cox模型分析

1.單變量Cox回歸

library("survival")
library("survminer")
res.cox <- coxph(Surv(time, status) ~ sex, data =  lung)
summary(res.cox)
res.cox <- coxph(Surv(time, status) ~ sex, data =  lung)
summary(res.cox)
summary的結(jié)果：

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

Cox回歸結(jié)果可以解釋如下：

1. 統(tǒng)計(jì)顯著性。標(biāo)記為“z”的列給出了Wald統(tǒng)計(jì)值。它對(duì)應(yīng)于每個(gè)回歸系數(shù)與其標(biāo)準(zhǔn)誤差的比率（z = coef / se（coef））。wald統(tǒng)計(jì)評(píng)估是否beta（ββ）系數(shù)在統(tǒng)計(jì)上顯著不同于0。從上面的輸出，我們可以得出結(jié)論，變量性別具有高度統(tǒng)計(jì)學(xué)意義的系數(shù)。

2. 回歸系數(shù)。Cox模型結(jié)果中要注意的第二個(gè)特征是回歸系數(shù)（coef）的符號(hào)。一個(gè)積極的信號(hào)意味著危險(xiǎn)（死亡風(fēng)險(xiǎn)）較高，因此對(duì)于那些變量值較高的受試者，預(yù)后更差。變量性被編碼為數(shù)字向量。1：男，2：女。Cox模型的R總結(jié)給出了第二組相對(duì)于第一組，即女性與男性的風(fēng)險(xiǎn)比（HR）。性別的β系數(shù)= -0.53表明在這些數(shù)據(jù)中，女性的死亡風(fēng)險(xiǎn)（低存活率）低于男性。

3. 危害比例。指數(shù)系數(shù)（exp（coef）= exp（-0.53）= 0.59）也稱為風(fēng)險(xiǎn)比，給出協(xié)變量的效應(yīng)大小。例如，女性（性別= 2）將危害降低了0.59倍，即41％。女性與預(yù)后良好相關(guān)。
   風(fēng)險(xiǎn)比的置信區(qū)間。總結(jié)結(jié)果還給出了風(fēng)險(xiǎn)比（exp（coef））的95％置信區(qū)間的上限和下限，下限95％界限= 0.4237，上限95％界限= 0.816。

4. 全球統(tǒng)計(jì)學(xué)意義的模型。最后，輸出為模型的總體顯著性提供了三個(gè)替代測(cè)試的p值：可能性比率測(cè)試，Wald測(cè)試和得分logrank統(tǒng)計(jì)。這三種方法是漸近等價(jià)的。對(duì)于足夠大的N，他們會(huì)得到相似的結(jié)果。對(duì)于小N來說，它們可能有所不同。似然比檢驗(yàn)對(duì)于小樣本量具有更好的表現(xiàn)，所以通常是優(yōu)選的。

2.多變量Cox回歸

要一次性將單變量coxph函數(shù)應(yīng)用于多個(gè)協(xié)變量，請(qǐng)輸入：

res.cox <- coxph(Surv(time, status) ~ age + sex + ph.ecog, data =  lung)
summary(res.cox)
Call:
coxph(formula = Surv(time, status) ~ age + sex + ph.ecog, data = lung)

  n= 227, number of events= 164
   (1 observation deleted due to missingness)

             coef exp(coef)  se(coef)      z Pr(>|z|)
age      0.011067  1.011128  0.009267  1.194 0.232416
sex     -0.552612  0.575445  0.167739 -3.294 0.000986 ***
ph.ecog  0.463728  1.589991  0.113577  4.083 4.45e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

        exp(coef) exp(-coef) lower .95 upper .95
age        1.0111     0.9890    0.9929    1.0297
sex        0.5754     1.7378    0.4142    0.7994
ph.ecog    1.5900     0.6289    1.2727    1.9864

Concordance= 0.637  (se = 0.025 )
Likelihood ratio test= 30.5  on 3 df,   p=1e-06
Wald test            = 29.93  on 3 df,   p=1e-06
Score (logrank) test = 30.5  on 3 df,   p=1e-06

所有3個(gè)整體測(cè)試（可能性，Wald 和 得分）的p值都是顯著的，表明該模型是顯著的。這些測(cè)試評(píng)估了所有的beta（ββ）為0.在上面的例子中，檢驗(yàn)統(tǒng)計(jì)是完全一致的，綜合零假設(shè)被完全拒絕。

在多變量Cox分析中，協(xié)變量性別和ph.ecog仍然顯著（p <0.05）。然而，協(xié)變量年齡并不顯著（P = 0.17，這比0.05）。

性別p值為0.000986，危險(xiǎn)比HR = exp（coef）= 0.58，表明患者性別和死亡風(fēng)險(xiǎn)降低之間有很強(qiáng)的關(guān)系。協(xié)變量的風(fēng)險(xiǎn)比可以解釋為對(duì)風(fēng)險(xiǎn)的倍增效應(yīng)。例如，保持其他協(xié)變量不變，女性（性別= 2）將危害降低0.58倍，即42％。我們的結(jié)論是，女性與良好的預(yù)后相關(guān)。

類似地，ph.ecog的p值是4.45e-05，危險(xiǎn)比HR = 1.59，表明ph.ecog值與死亡風(fēng)險(xiǎn)增加之間的強(qiáng)關(guān)系。保持其他協(xié)變量不變，ph.ecog值越高，生存率越差。

相比之下，年齡的p值現(xiàn)在是p = 0.23。風(fēng)險(xiǎn)比HR = exp（coef）= 1.01，95％置信區(qū)間為0.99至1.03。由于HR的置信區(qū)間為1，這些結(jié)果表明，年齡在調(diào)整了ph值和患者性別后對(duì)HR的差異的貢獻(xiàn)較小，并且僅趨向顯著性。例如，保持其他協(xié)變量不變，額外的年齡會(huì)導(dǎo)致每日死亡危險(xiǎn)因素為exp（beta）= 1.01或1％，這不是一個(gè)重要的貢獻(xiàn)。

可視化估計(jì)的生存時(shí)間分布

已經(jīng)將Cox模型擬合到數(shù)據(jù)中，可以在特定風(fēng)險(xiǎn)組的任何給定時(shí)間點(diǎn)可視化預(yù)測(cè)的存活比例。函數(shù)survfit（）估計(jì)生存比例，默認(rèn)情況下為協(xié)變量的平均值。

繪制生存曲線：

kmfit<-survfit(Surv(time, status) ~ sex, data =  lung)
ggsurvplot(kmfit,
              pval = TRUE, conf.int = TRUE,
              risk.table = TRUE, 
              risk.table.col = "strata", 
              linetype = "strata", 
              surv.median.line = "hv", 
              ggtheme = theme_bw(), 
              palette = c("#E7B800", "#2E9FDF"))

以上就是關(guān)于“R語言Cox與KM生存結(jié)果分析”這篇文章的內(nèi)容，相信大家都有了一定的了解，希望小編分享的內(nèi)容對(duì)大家有幫助，若想了解更多相關(guān)的知識(shí)內(nèi)容，請(qǐng)關(guān)注億速云行業(yè)資訊頻道。