生存分析基礎(chǔ)

生存分析對應(yīng)于一組統(tǒng)計(jì)方法，用于調(diào)查感興趣事件發(fā)生所花費(fèi)的時(shí)間。

生存分析可用于許多領(lǐng)域，例如：

• 用于患者生存時(shí)間分析的癌癥研究，
• “事件歷史分析”的社會(huì)學(xué)，
• 在工程中用于“故障時(shí)間分析”。

在癌癥研究中，典型的研究問題如下：

• 某些臨床特征對患者生存有何影響？
• 一個(gè)人生存3年的概率是多少？
• 兩組患者的生存率是否存在差異？

內(nèi)容

目標(biāo)

本章的目的是描述生存分析的基本概念。在癌癥研究中，大多數(shù)生存分析使用以下方法：

• Kaplan-Meier圖（Kaplan-Meier plots）可視化生存曲線
• 對數(shù)秩檢驗(yàn)（Log-rank test）以比較兩組或更多組的生存曲線
• 用Cox比例風(fēng)險(xiǎn)回歸（Cox proportional hazards regression）描述變量對生存的影響。下一章將討論Cox模型：Cox比例風(fēng)險(xiǎn)模型。

在這里，我們將從解釋生存分析的基本概念開始，包括：

• 如何生成和解釋生存曲線，
• 以及如何量化和測試兩組或更多組患者之間的生存差異。

然后，我們將繼續(xù)使用Cox比例風(fēng)險(xiǎn)模型描述多元分析。

基本概念

在這里，我們首先定義生存分析的基本術(shù)語，包括：

• 生存時(shí)間和事件（Survival time and event）
• 刪失（Censoring）
• 生存函數(shù)和風(fēng)險(xiǎn)函數(shù)（Survival function and hazard function）

癌癥研究中的生存時(shí)間和事件類型（Survival time and event）

有不同類型的事件，包括：

• 復(fù)發(fā)（Relapse）
• 進(jìn)展（Progression）
• 死亡（Death）

從“響應(yīng)到治療”（完全緩解）到發(fā)生關(guān)注事件的時(shí)間通常稱為生存時(shí)間（或事件發(fā)生時(shí)間）。癌癥研究中兩個(gè)最重要的標(biāo)度包括：i）死亡時(shí)間；ii）無復(fù)發(fā)生存時(shí)間，對應(yīng)于對治療的反應(yīng)與疾病復(fù)發(fā)之間的時(shí)間。也稱為無病生存時(shí)間（disease-free survival time）和無事件生存時(shí)間（event-free survival time）。

刪失（Censoring）

如上所述，生存分析著眼于直到發(fā)生感興趣事件（復(fù)發(fā)或死亡）之前的預(yù)期持續(xù)時(shí)間。但是，在研究期間內(nèi)某些人可能未觀察到該事件，從而產(chǎn)生了所謂的刪失（Censoring）現(xiàn)象。

審查可能以下列方式出現(xiàn)：

• 患者在研究時(shí)間段內(nèi)尚未（尚未）經(jīng)歷感興趣的事件，例如復(fù)發(fā)或死亡；
• 研究期間患者失去隨訪；
• 患者經(jīng)歷了另一種事件，因此無法進(jìn)行進(jìn)一步的隨訪。

在生存分析中會(huì)遇到這種現(xiàn)象，稱為右側(cè)刪失。

• 這里補(bǔ)充一下右側(cè)刪失和左側(cè)刪失的意思：以右側(cè)為例，當(dāng)患者發(fā)生上述情況時(shí)，在時(shí)間軸這個(gè)時(shí)間點(diǎn)的右側(cè)（即該時(shí)間點(diǎn)之后）數(shù)據(jù)點(diǎn)是缺失的，因此稱為右側(cè)刪失。臨床床上經(jīng)常遇到的是右側(cè)刪失。這樣左側(cè)刪失也容易理解了，既被隨訪者由于某些原因在時(shí)間軸內(nèi)某一點(diǎn)之前沒能參與隨訪，因此在改時(shí)間點(diǎn)之前（既時(shí)間軸左側(cè)）數(shù)據(jù)是缺失的，因此稱為左側(cè)刪失。
• 比如研究者想跟蹤調(diào)查青少年12歲至18歲之前的視力變化，如果某個(gè)被調(diào)查者在14歲才開始進(jìn)行隨訪就會(huì)產(chǎn)生左側(cè)缺失，如果某個(gè)被調(diào)查者在14歲由于玩游戲過度而住院無法繼續(xù)參與隨訪就會(huì)產(chǎn)生右側(cè)缺失。

生存和風(fēng)險(xiǎn)函數(shù)

使用兩個(gè)相關(guān)的概率來描述生存數(shù)據(jù)：生存概率和危險(xiǎn)概率。
生存函數(shù)，也被稱為幸存者函數(shù)S(t) 是從時(shí)間起源（例如癌癥診斷）到指定的未來時(shí)間t生存的概率。
風(fēng)險(xiǎn)函數(shù)，記為h(t)，是在時(shí)間t被觀察的個(gè)人發(fā)生某項(xiàng)事件的概率。
注意，生存函數(shù)側(cè)重于沒有事件，危害函數(shù)著重于事件發(fā)生。

Kaplan-Meier生存估計(jì)

Kaplan-Meier（KM）方法是一種非參數(shù)方法，用于根據(jù)觀察到的生存時(shí)間估算生存概率（Kaplan和Meier，1958年）。
時(shí)間點(diǎn) t(i) 的生存概率 計(jì)算如下：

估計(jì)的概率
估計(jì)概率（S（t））是僅在每個(gè)事件發(fā)生時(shí)改變值的階躍函數(shù)。也可以計(jì)算生存概率的置信區(qū)間。 KM生存曲線是KM生存概率與時(shí)間的關(guān)系曲線，它提供了一個(gè)有用的數(shù)據(jù)摘要，可用于估計(jì)中位生存時(shí)間等指標(biāo)。

R中的生存分析

安裝并加載所需的R包

我們將使用兩個(gè)R包：

• survival 用于計(jì)算存活分析
• survminer 用于總結(jié)和可視化生存分析結(jié)果
• 安裝軟件包

install.packages(c("survival", "survminer"))
library("survival")
library("survminer")
#我們將使用生存包中提供的肺癌數(shù)據(jù)。
data("lung")
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

• 機(jī)構(gòu)：機(jī)構(gòu)代碼
• 時(shí)間：以天為單位的生存時(shí)間
• 狀態(tài)：審查狀態(tài)1 =審查，2 =失效
• 年齡：歲
• 性別：男= 1女= 2
• ph.ecog：ECOG成績得分（0 =好5 =死）
• ph.karno：醫(yī)師對Karnofsky成績評分（差= 0-良好= 100）
• pat.karno：Karnofsky表現(xiàn)評分，按患者評分
• meal.cal：用餐時(shí)消耗的卡路里
• wt.loss：最近六個(gè)月的體重減輕

計(jì)算生存曲線：survfit（）

我們想按性別計(jì)算生存率。
生存軟件包中的功能survfit（）可用于計(jì)算kaplan-Meier生存估計(jì)。其主要論點(diǎn)包括：
使用Surv（）函數(shù)創(chuàng)建的生存對象

• 以及包含變量的數(shù)據(jù)集。
• 要計(jì)算生存曲線，請鍵入：

fit <- survfit(Surv(time, status) ~ sex, data = lung)
print(fit)

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

默認(rèn)情況下，函數(shù)print（）顯示生存曲線的簡短摘要。它打印觀察值，事件數(shù)，中位生存率和中位值的置信度限制。
如果要顯示生存曲線的更完整摘要，請鍵入以下內(nèi)容：

# Summary of survival curves
summary(fit)
# Access to the sort summary table
summary(fit)$table

訪問survfit（）返回的值

函數(shù)survfit（）返回變量列表，包括以下組件：

• n：每條曲線中的對象總數(shù)。
• 時(shí)間：曲線上的時(shí)間點(diǎn)。
• n。風(fēng)險(xiǎn)：時(shí)間t處有風(fēng)險(xiǎn)的受試者人數(shù)
• n.event：在時(shí)間t發(fā)生的事件數(shù)。
• n。審查者：在時(shí)間t退出事件而不發(fā)生風(fēng)險(xiǎn)的被審查者的數(shù)量。
• 下，上：曲線的置信度上限和下限。
• 分層：表示曲線估計(jì)的分層。如果地層不為NULL，則結(jié)果中有多條曲線。層次的水平（一個(gè)因子）是曲線的標(biāo)簽。
  可以按以下方式訪問組件

d <- data.frame(time = fit$time,
                  n.risk = fit$n.risk,
                  n.event = fit$n.event,
                  n.censor = fit$n.censor,
                  surv = fit$surv,
                  upper = fit$upper,
                  lower = fit$lower
                  )
head(d)

  time n.risk n.event n.censor      surv     upper     lower
1   11    138       3        0 0.9782609 1.0000000 0.9542301
2   12    135       1        0 0.9710145 0.9994124 0.9434235
3   13    134       2        0 0.9565217 0.9911586 0.9230952
4   15    132       1        0 0.9492754 0.9866017 0.9133612
5   26    131       1        0 0.9420290 0.9818365 0.9038355
6   30    130       1        0 0.9347826 0.9768989 0.8944820

可視化生存曲線

我們將使用功能ggsurvplot（）[在Survminer R包中]生成兩組受試者的生存曲線。
也可以顯示：

• 使用參數(shù)conf.int = TRUE對生存函數(shù)的95％置信度限制。
• 使用option risk.table按時(shí)間劃分處于風(fēng)險(xiǎn)中的個(gè)體的數(shù)量和/或百分比。risk.table的允許值包括：
  • RUE或FALSE指定是否顯示風(fēng)險(xiǎn)表。默認(rèn)值為FALSE。
  • “絕對”或“百分比”：分別顯示按時(shí)間劃分處于風(fēng)險(xiǎn)中的受試者的絕對數(shù)和百分比。使用“ abs_pct”顯示絕對數(shù)字和百分比。
• 對數(shù)秩檢驗(yàn)的p值，使用pval = TRUE比較組。
• 使用參數(shù)surv.median.line在中值生存時(shí)的水平/垂直線。允許的值包括c（“ none”，“ hv”，“ h”，“ v”）之一。v：垂直，h：水平。

# Change color, linetype by strata, risk.table color by strata
ggsurvplot(fit,
          pval = TRUE, conf.int = TRUE,
          risk.table = TRUE, # Add risk table
          risk.table.col = "strata", # Change risk table color by groups
          linetype = "strata", # Change line type by groups
          surv.median.line = "hv", # Specify median survival
          ggtheme = theme_bw(), # Change ggplot2 theme
          palette = c("#E7B800", "#2E9FDF"))

可以使用以下參數(shù)進(jìn)一步定制該圖：

• conf.int.style =“步驟”*以更改置信區(qū)間帶的樣式。
• xlab*更改x軸標(biāo)簽。
• break.time.by = 200*以時(shí)間間隔將x軸斷開200。
• risk.table =“ abs_pct”，*以顯示處于風(fēng)險(xiǎn)中的個(gè)人的絕對數(shù)量和百分比。
• risk.table.y.text.col = TRUE，risk.table.y.text = FALSE*在風(fēng)險(xiǎn)表圖例的文本注釋中提供條形而不是名稱。
• ncensor.plot = TRUE，*以繪制時(shí)間t處受檢對象的數(shù)量。正如Marcin Kosinski所建議的那樣，這是對生存曲線的一個(gè)很好的附加反饋，因此人們可以認(rèn)識到：生存曲線看起來如何，風(fēng)險(xiǎn)集的數(shù)量是多少，以及風(fēng)險(xiǎn)集變小的原因是什么？是由事件還是受審查事件引起的？
• legend.labs*更改圖例標(biāo)簽。

ggsurvplot(
   fit,                     # survfit object with calculated statistics.
   pval = TRUE,             # show p-value of log-rank test.
   conf.int = TRUE,         # show confidence intervals for 
                            # point estimaes of survival curves.
   conf.int.style = "step",  # customize style of confidence intervals
   xlab = "Time in days",   # customize X axis label.
   break.time.by = 200,     # break X axis in time intervals by 200.
   ggtheme = theme_light(), # customize plot and risk table with a theme.
   risk.table = "abs_pct",  # absolute number and percentage at risk.
  risk.table.y.text.col = T,# colour risk table text annotations.
  risk.table.y.text = FALSE,# show bars instead of names in text annotations
                            # in legend of risk table.
  ncensor.plot = TRUE,      # plot the number of censored subjects at time t
  surv.median.line = "hv",  # add the median survival pointer.
  legend.labs = 
    c("Male", "Female"),    # change legend labels.
  palette = 
    c("#E7B800", "#2E9FDF") # custom color palettes.
)

Kaplan-Meier圖可以解釋如下：

橫軸（x軸）表示以天為單位的時(shí)間，縱軸（y軸）表示生存的可能性或生存的人口比例。線代表兩組的存活曲線。曲線中的垂直下降表示事件。曲線上的垂直刻度線表示此時(shí)已對患者進(jìn)行檢查。

• 在零時(shí)，生存概率為1.0（或100％的參與者還活著）。
• 在時(shí)間250，性別= 1的存活概率約為0.55（或55％），性別= 2的存活概率約為- 0.75（或75％）。
• 性別= 2的中位生存期約為270天，性別= 2的中位生存期約為426天，這表明性別= 2的生存期高于性別= 1

可以使用以下代碼獲得每組的中位生存時(shí)間：

summary(fit)$table

      records n.max n.start events   *rmean *se(rmean) median 0.95LCL 0.95UCL
sex=1     138   138     138    112 325.0663   22.59845    270     212     310
sex=2      90    90      90     53 458.2757   33.78530    426     348     550

每組的中位生存時(shí)間代表生存概率S（t）為0.5的時(shí)間。
性別= 1（男性）的中位生存時(shí)間為270天，而性別= 2（女性）則為426天。與男性相比，女性肺癌似乎具有生存優(yōu)勢。但是，要評估這種差異是否具有統(tǒng)計(jì)學(xué)顯著性，需要進(jìn)行正式的統(tǒng)計(jì)檢驗(yàn)，這將在下一節(jié)中討論。

注意，置信極限在曲線的尾部很寬，因此很難進(jìn)行有意義的解釋。這可以通過以下事實(shí)來解釋：在實(shí)踐中，通常有一些患者在隨訪結(jié)束時(shí)迷失了隨訪或存活。因此，在隨訪結(jié)束之前在x軸上縮短圖可能是明智的（Pocock等，2002）。

可以使用參數(shù)xlim縮短生存曲線，如下所示：

ggsurvplot(fit,
          conf.int = TRUE,
          risk.table.col = "strata", # Change risk table color by groups
          ggtheme = theme_bw(), # Change ggplot2 theme
          palette = c("#E7B800", "#2E9FDF"),
          xlim = c(0, 600))

注意，可以使用參數(shù)fun指定三個(gè)經(jīng)常使用的轉(zhuǎn)換：

• “ log”：幸存者功能的對數(shù)轉(zhuǎn)換，
• “事件”：繪制累積事件（f（y）= 1-y）。也稱為累積發(fā)生率
• “ cumhaz”繪制了累積危害函數(shù)（f（y）= -log（y））

例如，要繪制累積事件，請鍵入：

ggsurvplot(fit,
          conf.int = TRUE,
          risk.table.col = "strata", # Change risk table color by groups
          ggtheme = theme_bw(), # Change ggplot2 theme
          palette = c("#E7B800", "#2E9FDF"),
          fun = "event")

累積性危險(xiǎn)是常用來估計(jì)危險(xiǎn)概率。定義為：

H(t)=?log(survivalfunction)=?log(S(t))

累積危險(xiǎn)（H（t））可以解釋為死亡的累積力。換言之，如果事件是可重復(fù)的過程，則它對應(yīng)于時(shí)間t時(shí)每個(gè)個(gè)體預(yù)期的事件數(shù)。

要繪制累積危害，請鍵入以下內(nèi)容：

ggsurvplot(fit,
          conf.int = TRUE,
          risk.table.col = "strata", # Change risk table color by groups
          ggtheme = theme_bw(), # Change ggplot2 theme
          palette = c("#E7B800", "#2E9FDF"),
          fun = "cumhaz")

Kaplan-Meier生命表：生存曲線摘要

如上所述，您可以使用函數(shù)summary（）來獲得生存曲線的完整摘要：

summary(fit)

還可以使用功能surv_summary（）[在survminer程序包中]獲取生存曲線的摘要。與默認(rèn)的summary（）函數(shù)相比，surv_summary（）創(chuàng)建一個(gè)數(shù)據(jù)框，其中包含來自survfit結(jié)果的不錯(cuò)的摘要。

res.sum <- surv_summary(fit)
head(res.sum)
  time n.risk n.event n.censor      surv    std.err     upper     lower strata sex
1   11    138       3        0 0.9782609 0.01268978 1.0000000 0.9542301  sex=1   1
2   12    135       1        0 0.9710145 0.01470747 0.9994124 0.9434235  sex=1   1
3   13    134       2        0 0.9565217 0.01814885 0.9911586 0.9230952  sex=1   1
4   15    132       1        0 0.9492754 0.01967768 0.9866017 0.9133612  sex=1   1
5   26    131       1        0 0.9420290 0.02111708 0.9818365 0.9038355  sex=1   1
6   30    130       1        0 0.9347826 0.02248469 0.9768989 0.8944820  sex=1   1

函數(shù)surv_summary（）返回包含以下列的數(shù)據(jù)幀：

• 時(shí)間：曲線有階躍的時(shí)間點(diǎn)。
• n。風(fēng)險(xiǎn)：處于t風(fēng)險(xiǎn)的受試者人數(shù)。
• n.event：在時(shí)間t發(fā)生的事件數(shù)。
• n.censor：審查事件的數(shù)量。
• surv：估計(jì)生存概率。
• std.err：生存標(biāo)準(zhǔn)誤差。
• upper：置信區(qū)間的上限
• 較低：置信區(qū)間的下限
• 分層：表示曲線估計(jì)的分層。層次的水平（一個(gè)因子）是曲線的標(biāo)簽。

在生存曲線已擬合一個(gè)或多個(gè)變量的情況下，surv_summary對象包含表示變量的額外列。這使得可以按層次或某些因素組合來考慮ggsurvplot的輸出。

surv_summary對象還有一個(gè)名為“表”的屬性，其中包含有關(guān)生存曲線的信息，包括具有置信區(qū)間的生存中位數(shù)，以及受試者總數(shù)和每條曲線中的事件數(shù)。要訪問屬性“表”，請輸入以下命令：

attr(res.sum, "table")

比較生存曲線的Log-Rank檢驗(yàn)：survdiff（）

log-rank test 是比較兩條或更多條生存曲線的最廣泛使用的方法。零假設(shè)是兩組之間的生存率沒有差異。對數(shù)秩檢驗(yàn)是一種非參數(shù)檢驗(yàn)，它不對生存分布做出任何假設(shè)。本質(zhì)上，對數(shù)秩檢驗(yàn)將每個(gè)組中觀察到的事件數(shù)與如果原假設(shè)為真（即，生存曲線相同）時(shí)的預(yù)期事件數(shù)進(jìn)行比較。對數(shù)秩統(tǒng)計(jì)量近似作為卡方檢驗(yàn)統(tǒng)計(jì)量分布。

函數(shù)survdiff（）[在生存包中]可用于比較兩個(gè)或更多生存曲線的對數(shù)秩檢驗(yàn)。

survdiff（）可以如下使用：

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/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 

該函數(shù)返回組件列表，包括：

• n：每組的科目數(shù)。
• obs：每組中事件的加權(quán)觀測數(shù)量。
• exp：每個(gè)組中加權(quán)的預(yù)期事件數(shù)。
• chisq：用于檢驗(yàn)相等性的卡方統(tǒng)計(jì)量。
• 階層：（可選）每個(gè)階層中包含的主題數(shù)。

生存差異的對數(shù)秩檢驗(yàn)給出p值為p = 0.0013，表明性別群體的生存差異顯著。

擬合復(fù)雜的生存曲線

在本節(jié)中，我們將使用多個(gè)因素的組合來計(jì)算生存曲線。接下來，我們將結(jié)合多種因素來研究ggsurvplot（）的輸出

1. 使用結(jié)腸數(shù)據(jù)集擬合（復(fù)雜）生存曲線

require("survival")
fit2 <- survfit( Surv(time, status) ~ sex + rx + adhere,
                data = colon )

2. 使用survminer可視化輸出。下圖顯示了根據(jù)rx＆粘附力值按性別分面的生存曲線。

# Plot survival curves by sex and facet by rx and adhere
ggsurv <- ggsurvplot(fit2, fun = "event", conf.int = TRUE,
                     ggtheme = theme_bw())

ggsurv$plot +theme_bw() + 
  theme (legend.position = "right")+
  facet_grid(rx ~ adhere)

概要

生存分析是用于數(shù)據(jù)分析的一組統(tǒng)計(jì)方法，其中感興趣的結(jié)果變量是事件發(fā)生之前的時(shí)間。
生存數(shù)據(jù)通常根據(jù)兩個(gè)相關(guān)功能進(jìn)行描述和建模：

• 幸存者函數(shù)代表一個(gè)人從起源時(shí)間到超過時(shí)間t的某個(gè)時(shí)間生存的概率。通常用Kaplan-Meier方法估算。對數(shù)秩檢驗(yàn)可用于測試組（例如治療組）的生存曲線之間的差異。

• 危害函數(shù)給出了一次事件的瞬時(shí)可能性，并給出了到那個(gè)時(shí)間的生存率。它主要用作診斷工具或用于指定生存分析的數(shù)學(xué)模型。

在本文中，我們演示了如何使用兩個(gè)R包（survival（用于分析）和 survminer（用于可視化））的組合來執(zhí)行和可視化生存分析。

轉(zhuǎn)自：http://www.itdecent.cn/p/9838e8004c5d