1. 程式人生 > 其它 >【視訊】R語言生存分析原理與晚期肺癌患者分析案例|資料分享

【視訊】R語言生存分析原理與晚期肺癌患者分析案例|資料分享

原文連結:http://tecdat.cn/?p=10278

原文出處:拓端資料部落公眾號

生存分析(也稱為工程中的可靠性分析)的目標是在協變數和事件時間之間建立聯絡。生存分析的名稱源於臨床研究,其中預測死亡時間,即生存,通常是主要目標。

視訊:R語言生存分析原理與晚期肺癌患者分析案例

R語言生存分析Survival analysis原理與晚期肺癌患者分析案例

,時長08:41

生存分析是一種迴歸問題(人們想要預測一個連續值),但有一個轉折點。它與傳統迴歸的不同之處在於,在生存分析中,結果變數既有一個事件,也有一個與之相關的時間值,部分訓練資料只能被部分觀察——它們是被刪失的。本文用R語言生存分析晚期肺癌患者資料

檢視文末了解資料獲取方式)。

普通最小二乘迴歸方法不足,因為事件發生的時間通常不是正態分佈的,並且模型無法處理刪失,但這在生存資料中很常見。

為什麼要做生存分析:右刪失

在某些情況下,可能無法觀察到事件時間:這通常稱為 右刪失。在以死亡為事件的臨床試驗中,當發生以下情況之一時,就會發生這種情況。1。當一定數量的參與者死亡時,研究結束。2。參與者退出研究。3。 研究達到預定的結束時間,並且一些參與者存活到結束。在每種情況下,倖存的參與者離開研究後,我們都不知道他們會發生什麼。然後我們有一個問題:

當對於某些個體,我們只觀察到他們的事件時間的下限時,我們如何對經驗分佈進行建模或進行非負迴歸?

上圖說明了右刪失。對於參與者 1,我們看到他們何時死亡。參與者 2 退出了,我們知道他們一直活到那時,但不知道後來發生了什麼。對於參與者 3,我們知道他們活到了預定的研究結束,但又不知道之後發生了什麼。

生存函式和風險函式

生存分析中的兩個關鍵工具是生存函式和風險函式。

生存函式:它是一個函式,用於給出我們有興趣知道的任何物件是否會在任何指定時間之後存活的概率。在數學上它可以由以下公式表示 

其中 S(t) 是一個生存函式,其中 T 是一個連續隨機變數,是一個事件的時間。F(t) 是區間[0,∞) 上的累積分佈函式。

我們也可以用風險函式來寫生存函式。假設事件尚未發生 ,風險率λ(t) 是事件在時間t發生的瞬時概率的主要值。

那麼關鍵問題是如何估計風險和/或生存函式。

Kaplan Meier的非引數估計

在非引數生存分析中,我們要估計生存函式沒有協變數,並且有刪失。如果我們沒有刪失,我們可以從經驗 CDF 開始. 這個等式簡潔地表示:

有多少人隨著時間的推移而死亡? 那麼生存函式就是:還有多少人還活著?

但是,我們無法回答一些人被時間t刪失時提出的這個問題.

雖然我們不一定知道有多少人在任意時間t倖存下來,我們知道研究中有多少人仍然處於風險之中。我們可以使用它來代替。將學習時間劃分區間, 其中每個ti是參與者的事件時間或刪失時間。假設參與者只能在觀察到的事件時間失效。假設沒有人在同一時間死去(沒有關係),我們可以檢視每次有人死去的時間。我們說在那個特定時間死亡的概率是,並說在任何其他時間死亡的概率是0.

在溫和的假設下,包括參與者具有獨立且相同分佈的事件時間,並且刪失和事件時間是獨立的,這給出了一個一致的估計量。上圖給出了一個簡單案例的 Kaplan Meier 估計示例。

生存分析用於各種領域

例如:

  • 用於患者生存時間分析的癌症研究,

  • “事件歷史分析”的社會學,

  • 在工程中用於“故障時間分析”。

在癌症研究中,典型的研究問題如下:

  • 某些臨床特徵對患者生存有何影響

  • 一個人能活3年的概率是多少?

  • 患者組之間的生存率是否存在差異?

第1部分:生存分析簡介 

本簡報將介紹生存分析 ,參考:

Clark, T., Bradburn, M., Love, S., & Altman, D. (2003). Survival analysis part I: Basic concepts and first analyses. 232-238. ISSN 0007-0920.

我們今天將使用的一些軟體包包括:

  • lubridate

  1.   library(survival)
  2.    

什麼是生存資料?

事件時間資料由不同的開始時間和結束時間組成。

癌症的例子

  • 從手術到死亡的時間

  • 從治療開始到進展的時間

  • 從響應到復發的時間

其他領域的例子

事件發生時間資料在許多領域都很常見,包括但不限於

  • 從艾滋病毒感染到艾滋病發展的時間

  • 心臟病發作的時間

  • 藥物濫用發生的時間

  • 機器故障時間

生存分析別名

由於生存分析在許多其他領域很常見,因此也有其他名稱

  • 可靠性分析

  • 持續時間分析

  • 事件歷史分析

  • 事件發生時間分析

肺資料集

資料包含來自北中部癌症治療組的晚期肺癌患者。今天我們將用來演示方法的一些變數包括

  • 時間:以天為單位的生存時間

  • 狀態:刪失狀態1 =刪失,2 =失效

  • 性別:男= 1女= 2

刪失型別

某個主題可能由於以下原因而被刪失:

  • 後續損失

  • 退出研究

  • 固定學習期結束前沒有活動

具體來說,這些是刪失的示例。

分配隨訪時間

  • 受刪失的主題仍會提供資訊,因此必須適當地包含在分析中

  • 隨訪時間的分佈存在偏差,在接受檢查的患者和有事件的患者之間可能有所不同

生存資料的組成部分

對於主題ii:

  • 活動時間Ti

  • 刪失時間Ci

  • 事件指標δi:

  • 1,如果觀察到的事件(即  Ti≤CiTi≤Ci)

  • 如果檢查,則為0(即  Ti>CiTi>Ci)

  1. 觀測時間Yi=min(Ti,Ci)Yi=min(Ti,Ci)

lung資料中提供了觀察時間和事件指示

  • 時間:以天為單位的生存時間(YiYi)

  • 狀態:刪失狀態1 =刪失,2 =死亡(δiδi)

在R中處理日期

資料通常帶有開始日期和結束日期,而不是預先計算的生存時間。第一步是確保將這些格式設定為R中的日期。

讓我們建立一個小的示例資料集,其中sx_date包含手術日期和last_fup_date上次隨訪日期的變數。

  1.   date_ex <- 
  2.     tibble(
  3.       sx_date = c("2007-06-22", "2004-02-13", "2010-10-27"), 
  4.       last\_fup\_date = c("2017-04-15", "2018-07-04", "2016-10-31")
  5.       )
  6.    
  7.   date_ex
  1.   ## # A tibble: 3 x 2
  2.   ##   sx\_date    last\_fup_date
  3.   ##   <chr>      <chr>        
  4.   ## 1 2007-06-22 2017-04-15   
  5.   ## 2 2004-02-13 2018-07-04   
  6.   ## 3 2010-10-27 2016-10-31

我們看到它們都是字元變數,通常都是這種情況,但是我們需要將它們格式化為日期。

格式化日期-基數R

  1.   date_ex %>% 
  2.     mutate(
  3.       sx\_date = as.Date(sx\_date, format = "%Y-%m-%d"), 
  4.       last\_fup\_date = as.Date(last\_fup\_date, format = "%Y-%m-%d") 
  5.       )
  1.   ## # A tibble: 3 x 2
  2.   ##   sx\_date    last\_fup_date
  3.   ##   <date>     <date>       
  4.   ## 1 2007-06-22 2017-04-15   
  5.   ## 2 2004-02-13 2018-07-04   
  6.   ## 3 2010-10-27 2016-10-31
  • 請注意,R格式必須包含分隔符和符號。例如,如果您的日期格式為m / d / Y,則需要format = "%m/%d/%Y"

格式化日期-lubridate程式包

我們還可以使用該lubridate包來格式化日期。在這種情況下,請使用ymd功能

  1.   date_ex %>% 
  2.     mutate(
  3.       sx\_date = ymd(sx\_date), 
  4.       last\_fup\_date = ymd(last\_fup\_date)
  5.       )
  1.   ## # A tibble: 3 x 2
  2.   ##   sx\_date    last\_fup_date
  3.   ##   <date>     <date>       
  4.   ## 1 2007-06-22 2017-04-15   
  5.   ## 2 2004-02-13 2018-07-04   
  6.   ## 3 2010-10-27 2016-10-31
  • 請注意,與基本R選項不同,不需要指定分隔符

計算生存時間 

現在日期已格式化,我們需要以某些單位(通常是幾個月或幾年)計算開始時間和結束時間之間的差。在base中R,用於difftime計算兩個日期之間的天數,然後使用將其轉換為數字值as.numeric。然後將除以365.25年的平均天數轉換為年。

  1.   date_ex %>% 
  2.     mutate(
  3.       os_yrs = 
  4.         as.numeric(
  5.           difftime(last\_fup\_date, 
  6.                    sx_date, 
  7.                    units = "days")) / 365.25
  8.       )
  1.   ## # A tibble: 3 x 3
  2.   ##   sx\_date    last\_fup\_date os\_yrs
  3.   ##   <date>     <date>         <dbl>
  4.   ## 1 2007-06-22 2017-04-15      9.82
  5.   ## 2 2004-02-13 2018-07-04     14.4 
  6.   ## 3 2010-10-27 2016-10-31      6.01

計算生存時間 

操作員可以%--%指定一個時間間隔,然後使用將該時間間隔轉換為經過的秒數as.duration,最後除以dyears(1),將其轉換為年數,從而得出一年中的秒數。

  1.   ## # A tibble: 3 x 3
  2.   ##   sx\_date    last\_fup\_date os\_yrs
  3.   ##   <date>     <date>         <dbl>
  4.   ## 1 2007-06-22 2017-04-15      9.82
  5.   ## 2 2004-02-13 2018-07-04     14.4 
  6.   ## 3 2010-10-27 2016-10-31      6.02

事件標標

對於生存資料的組成部分,我提到了事件指示器:

事件指標δiδi:

  • 1,如果觀察到的事件(即  Ti≤CiTi≤Ci)

  • 如果檢查,則為0(即  Ti>CiTi>Ci)

lung資料中,我們有:

  • 狀態:刪失狀態1 =刪失,2 =失效

生存函式

受試者可以存活超過指定時間的概率

S(t)=Pr(T>t)=1−F(t)S(t)=Pr(T>t)=1−F(t)

S(t)S(t):生存函式F(t)=Pr(T≤t)F(t)=Pr(T≤t):累積分佈函式

理論上,生存函式是平滑的;在實踐中,我們以離散的時間尺度觀察事件。

生存概率

  • 生存概率在某個時間,S(t)S(t),是存活超過該時間,考慮到個體已存活剛剛在此之前,時間的條件概率。

  • 可以估計為當時活著但沒有損失的隨訪患者人數除以當時的活著患者人數

  • 生存概率的Kaplan-Meier估計是這些條件概率的乘積

  • 在時間0,生存概率為1,即  S(t0)=1S(t0)=1

建立生存物件

Kaplan-Meier方法是估計生存時間和概率的最常用方法。這是一種非引數方法,可產生階躍函式,每次事件發生時,階躍下降。

  • 建立一個生存物件。對於每個主題,將有一個條目作為生存時間,+如果主題是經過刪失的,則後面跟一個。讓我們看一下前10個觀察值:

##  \[1\]  306   455  1010+  210   883  1022+  310   361   218   166

用Kaplan-Meier方法估算生存曲線

  • survfit函式根據公式建立生存曲線。讓我們為整個同類群組生成總體生存曲線,將其分配給object f1,然後檢視names該物件的:

names(f1)
  1.   ##  \[1\] "n"          "time"       "n.risk"     "n.event"    "n.censor"  
  2.   ##  \[6\] "surv"       "std.err"    "cumhaz"     "std.chaz"   "start.time"
  3.   ## \[11\] "type"       "logse"      "conf.int"   "conf.type"  "lower"     
  4.   ## \[16\] "upper"      "call"

survfit物件將用於建立生存曲線的一些關鍵元件包括:

  • time,其中包含每個時間間隔的起點和終點

  • surv,其中包含每個對應的生存概率 time

Kaplan-Meier圖 

現在, 繪製物件 獲得Kaplan-Meier圖。

  1.   plot(survfit(Surv(time, status) ~ 1, data = lung), 
  2.    
  • 基數R中的預設圖顯示了具有相關置信區間(虛線)的階躍函式(實線)

  • 水平線代表間隔的生存時間

  • 時間間隔由事件終止

  • 垂直線的高度顯示累積概率的變化

  • 帶有刻度線的經過刪失的觀察結果會減少間隔之間的累積生存期。

Kaplan-Meier圖 

建立在上ggplot2,並可用於建立Kaplan-Meier圖。

  • 預設圖  帶相關置信帶(陰影區域)的階躍函式(實線)。

  • 預設情況下,顯示了被檢查患者的刻度線,在此示例中,該刻度線本身有些模糊,可以使用選項將其取消 censor = FALSE

估計xx年生存

生存分析中經常需要關注的一個數量是生存超過一定數量(xx)年的概率。

例如,要估算生存到11年的可能性

  1.   ## Call: survfit(formula = Surv(time, status) ~ 1, data = lung)
  2.   ## 
  3.   ##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
  4.   ##   365     65     121    0.409  0.0358        0.345        0.486

我們發現本研究中11年生存的機率是41%。

同時顯示95%置信區間的相關上下限。

xx年生存率和生存曲線

11年存活率概率為在y軸上的點對應於11一年x軸的生存曲線。

Xx年生存率常常被錯誤估計

如果 使用“天真”的估計會怎樣?

228名患者中的121名到1年時死亡,因此:

-當 忽略42名患者在1年之前受到檢查的事實時, 會錯誤估計1個1個年生存率。

  • 正確的估計生存概率-年為41%。

忽略刪失對xx年生存的影響

  • 想象兩個研究,每個研究228個主題。每個研究中有165人死亡。一個沒有檢查(紅色線),63個病人被另一個(藍色線)檢查

  • 忽略刪失會導致總體生存概率被高估,因為被刪失的受試者僅在部分隨訪時間內提供資訊,然後落入風險範圍之外,從而降低了生存的累積概率

估計中位生存時間

生存分析中經常需要關注的另一個數量是平均生存時間,我們使用中位數對其進行量化。預計生存時間不會呈正態分佈,因此平均值不是適當的總結。

  1.   ## Call: survfit(formula = Surv(time, status) ~ 1, data = lung)
  2.   ## 
  3.   ##       n  events  median 0.95LCL 0.95UCL 
  4.   ##     228     165     310     285     363

我們看到中位生存時間為310天。還會顯示95%置信區間的上限和下限。

中位生存時間和生存曲線

中位生存時間是生存概率為0.50 

中位生存率常常被錯誤估計

總結165例死亡患者的中位生存時間

  1.   ##   median_surv
  2.   ## 1         226
  • 當 忽略被檢查患者也有助於隨訪的事實時, 會得出錯誤的估計中值生存時間226天。

  • 正確的中位生存時間估計是310天。

忽略刪失對中位數生存率的影響

  • 忽略刪失會造成人為降低的生存曲線,因為排除了受刪失患者貢獻的隨訪時間(紫色線)

  • 資料的真實生存曲線以lung藍色顯示,以進行比較

比較各組之間的生存時間

  • 我們可以使用對數秩檢驗進行組間重要性檢驗

  • 對數秩檢驗在整個隨訪時間內平均權衡觀察結果,是比較組間生存時間的最常用方法

  • 根據研究問題,有些版本可能會更重視早期或後期的隨訪,可能更合適

我們使用 函式獲得對數秩p值。例如,我們可以根據lung資料中的性別測試是否存在生存時間差異

  1.   ## Call:
  2.   ## 
  3.   ##         N Observed Expected (O-E)^2/E (O-E)^2/V
  4.   ## sex=1 138      112     91.6      4.55      10.3
  5.   ## sex=2  90       53     73.4      5.68      10.3
  6.   ## 
  7.   ##  Chisq= 10.3  on 1 degrees of freedom, p= 0.001

從survdiff物件中提取資訊

從 結果中提取p值

1 - pchisq(sd$chisq, length(sd$n) - 1)
## \[1\] 0.001311165

返回格式化的p值

## \[1\] 0.001

Cox迴歸模型

我們可能想量化單個變數的效應大小,或者將多個變數包括在迴歸模型中以說明多個變數的效應。

Cox迴歸模型是半引數模型,可用於擬合具有生存結果的單變數和多變量回歸模型。

h(t)h(t):危險或事件發生的瞬時速率h0(t)h0(t):基本基準危險

該模型的一些關鍵假設:

  • 非資訊刪失

  • 比例危險

_注意_:也可以使用用於生存結果的引數迴歸模型,但是本培訓將不涉及這些模型。

我們可以使用coxph函式擬合生存資料的迴歸模型,該函式Surv在左側使用一個物件,而在右側具有用於迴歸公式的標準語法R

  1.   ## Call:
  2.   ## 
  3.   ##        coef exp(coef) se(coef)      z       p
  4.   ## sex -0.5310    0.5880   0.1672 -3.176 0.00149
  5.   ## 
  6.   ## Likelihood ratio test=10.63  on 1 df, p=0.001111
  7.   ## n= 228, number of events= 165

格式化Cox迴歸結果

我們可以看到輸出的整潔版本broom

或使用

危險比

  • 來自Cox迴歸模型的關注數量是危險比(HR)。HR表示在任何特定時間點兩組之間的危險比率。

  • HR被解釋為感興趣事件中那些仍處於事件風險中的事件的瞬時發生率。

  • 如果您有一個迴歸引數ββ(來自estimate我們的列coxph),則HR = 經驗值(β)經驗值⁡(β)。

  • HR <1表示死亡危險降低,而HR> 1表示死亡危險增加。

  • 因此,我們的HR = 0.59意味著在任何給定時間,女性死亡的人數大約是男性的0.6倍。

第2部分:地標分析和時間相關協變數

在第1部分中,我們介紹了使用對數秩檢驗和Cox迴歸來檢驗感興趣的協變數與生存結果之間的關聯。

示例:腫瘤反應

示例:從治療開始就測量總生存期,關注的是對治療的完全反應與生存之間的關聯。

  • Anderson等人(JCO,1983)描述了在這種情況下,為什麼傳統方法(如對數秩檢驗或Cox迴歸)偏向於響應者,並提出了劃時代的方法。

  • 界標方法中的零假設是,從界標生存的過程不依賴於界標的響應狀態。

Anderson, J., Cain, K., & Gelber, R. (1983). Analysis of survival by tumor response. Journal of Clinical Oncology : Official Journal of the American Society of Clinical Oncology, 1(11), 710-9.

其他例子

癌症研究中可能尚未關注的其他一些可能的協變數包括:

  • 移植失敗

  • 移植物抗宿主病

  • 第二次切除

  • 輔助治療

  • 合規

  • 不良事件

示例資料 

137例骨髓移植患者的資料。變數包括:

  • T1 死亡時間或最後一次隨訪時間(天)

  • delta1死亡指標;1死0活

  • TA 急性移植物抗宿主病的時間(以天為單位)

  • deltaA急性移植物抗宿主病指標;1-發展為急性移植物抗宿主病,0-從未發展為急性移植物抗宿主病

讓我們載入資料以供整個示例使用

地標法

  1. 選擇基線之後的固定時間作為界標時間。_注意_:應在檢查資料之前根據臨床資訊進行操作

  2. 那些人群的子集至少跟蹤到里程碑時間。_注意_:請務必在地標時間之前報告由於關注或刪失事件而排除的號碼。

  3. 計算具有里程碑意義的時間,並應用傳統的對數秩檢驗或Cox迴歸

BMT資料感興趣的是急性移植物抗宿主病(aGVHD)和存活之間的關聯。但是aGVHD是在移植後進行評估的,這是我們的基線,也就是後續隨訪的開始時間。

步驟1選擇地標時間

通常,aGVHD發生在移植後的前90天內,因此我們使用90天的界標。

人們對急性移植物抗宿主病(aGVHD)與生存之間的關係感興趣。但是aGVHD是在移植後進行評估的,這是我們的基線,也就是後續隨訪的開始時間。

第2步:至少跟蹤到里程碑時間之前的人群的子集

這將我們的樣本量從137減少到122。

  • 所有15位被排除的患者均在90天里程碑之前死亡

人們對急性移植物抗宿主病(aGVHD)與生存之間的關係感興趣。但是aGVHD是在移植後進行評估的,這是我們的基線,也就是後續隨訪的開始時間。

步驟3根據地標計算隨訪時間,並應用傳統方法。

使用BMT資料的Cox迴歸界標示例

在Cox迴歸中, 可以使用中的subset選項coxph來排除那些在標誌性時間內沒有被隨訪的患者

時間相關協變數

界標分析的替代方法是合併時間相關的協變數。這可能更適合

  • 協變數的值隨時間變化

  • 沒有明顯的里程碑時間

時間相關協變數資料設定

對時間相關協變數的分析R需要建立特殊的資料集。 

BMT資料中沒有ID變數,這是建立特殊資料集所必需的,因此請建立一個名為的變數my_id

tmerge函式與event和函式一起使用tdc可建立特殊資料集。

  • tmerge 為每個患者的不同協變數值建立一個具有多個時間間隔的長資料集

  • event 建立新的事件指示器,以與新建立的時間間隔一致

  • tdc 建立與時間相關的協變數指標,以與新建立的時間間隔一致

時間相關協變數-單例患者

要了解其作用,讓我們看一下前5名患者的資料。

  1.    ##   my_id   T1 delta1   TA deltaA
  2.   ## 1     1 2081      0   67      1
  3.   ## 2     2 1602      0 1602      0
  4.   ## 3     3 1496      0 1496      0
  5.   ## 4     4 1462      0   70      1
  6.   ## 5     5 1433      0 1433      0

這些相同患者的新資料集

  1.   ##   my_id   T1 delta1 id tstart tstop death agvhd
  2.   ## 1     1 2081      0  1      0    67     0     0
  3.   ## 2     1 2081      0  1     67  2081     0     1
  4.   ## 3     2 1602      0  2      0  1602     0     0
  5.   ## 4     3 1496      0  3      0  1496     0     0
  6.   ## 5     4 1462      0  4      0    70     0     0
  7.   ## 6     4 1462      0  4     70  1462     0     1
  8.   ## 7     5 1433      0  5      0  1433     0     0

時間相關協變數-Cox迴歸

現在,我們可以分析這個時間依賴性協照常使用Cox迴歸與coxph 

摘要

我們發現,使用標誌性分析或時間依賴性協變數,急性移植物抗宿主病與死亡無顯著相關性。

通常,人們會希望使用地標分析對單個協變數進行視覺化, 使用帶有時間相關協變數的Cox迴歸進行單變數和多變數建模。

第3部分:競爭風險

什麼是競爭風險?

當物件在事件發生時間設定中發生多個可能的事件時

例子:

  • 復發

  • 因疾病死亡

  • 因其他原因死亡

  • 治療反應

在任何給定的研究中,所有這些(或其中一些 以及其他)可能都是可能的事件。

所以有什麼問題?

事件時間之間未觀察到的依賴性是導致需要特殊考慮的基本問題。

例如,可以想象復發的患者更有可能死亡,因此復發時間和死亡時間將不是獨立事件。

競爭風險的背景

存在多種潛在結果時的兩種分析方法:

  1. 給定事件的特定於原因的危險:這表示未因其他事件而失敗的事件中事件的每單位時間的發生率

  2. 給定事件的累積發生率:這表示事件每單位時間的發生率以及競爭事件的影響

這些方法中的每一種都可能僅闡明資料的一個重要方面,而有可能使其他方面難以理解,因此所選的方法應取決於感興趣的問題。

 黑色素瘤資料示例

它包含變數:

  • time 生存時間以天為單位,可能經過刪失。

  • status 1例死於黑色素瘤,2例存活,3例因其他原因死亡。

  • sex 1 =男性,0 =女性。

  • age 年歲。

  • year 操作。

  • thickness 腫瘤厚度(毫米)。

  • ulcer 1 =存在,0 =不存在。

黑色素瘤資料的累積發生率

在競爭風險的背景下估算累積發生率。

  1.   ## Estimates and Variances:
  2.   ## $est
  3.   ##           1000       2000       3000      4000      5000
  4.   ## 1 1 0.12745714 0.23013963 0.30962017 0.3387175 0.3387175
  5.   ## 1 3 0.03426709 0.05045644 0.05811143 0.1059471 0.1059471
  6.   ## 
  7.   ## $var
  8.   ##             1000         2000         3000        4000        5000
  9.   ## 1 1 0.0005481186 0.0009001172 0.0013789328 0.001690760 0.001690760
  10.   ## 1 3 0.0001628354 0.0002451319 0.0002998642 0.001040155 0.001040155

繪製累積發生率-基數R

生成 預設值的基本圖。

plot(ci_fit)

繪製累積發生率 

比較組之間的累積發生率

用於組間測試。

例如,Melanoma根據ulcer潰瘍的存在與否比較結果。測試結果可以在中找到Tests

ci_ulcer\[\["Tests"\]\]
  1.   ##        stat           pv df
  2.   ## 1 26.120719 3.207240e-07  1
  3.   ## 3  0.158662 6.903913e-01  1

按組繪製累積發生率 

按組繪製累積發生率-手動

_請注意,_我個人發現該ggcompetingrisks功能缺少自定義功能,尤其是與相比ggsurvplot。我通常會自己做圖,首先建立cuminc擬合結果的整潔資料集,然後再繪製結果。有關底層程式碼的詳細資訊,請參見此簡報的

繪製單個事件型別 

通常,只有一種型別的事件會引起人們的興趣,儘管我們仍要考慮競爭事件。在那種情況下,感興趣的事件可以單獨繪製。同樣,我首先通過建立cuminc擬合結果的整潔資料集,然後繪製結果來手動執行此操作。有關底層程式碼的詳細資訊,請參見此簡報的原始碼。

在風險表中新增數字

您可能想將風險表的數量新增到累積發生率圖中,而據我所知,沒有簡單的方法可以做到這一點。請參閱此簡報的原始碼中的一個示例

競爭風險迴歸

兩種方法:

  1. 特定原因風險

  • 當前沒有事件的受試者中給定事件型別的瞬時發生率

  • 使用Cox迴歸估算

  1. Subdistribution子分佈風險

  • 給定型別事件在沒有經歷過此類事件的受試者中的瞬時發生率

  • 使用Fine-Gray迴歸估算

黑色素瘤資料中的競爭風險迴歸-子分佈風險法Subdistribution

假設我們有興趣研究年齡和性別對黑色素瘤死亡的影響,而其他原因的死亡則是競爭事件。

  • crr 需要指定協變數作為矩陣

  • 如果 關注多個事件,則可以使用failcode選項請求其他事件的結果,預設情況下會返回failcode = 1

shr_fit
  1.   ## convergence:  TRUE 
  2.   ## coefficients:
  3.   ##     sex     age 
  4.   ## 0.58840 0.01259 
  5.   ## standard errors:
  6.   ## \[1\] 0.271800 0.009301
  7.   ## two-sided p-values:
  8.   ##  sex  age 
  9.   ## 0.03 0.18

在上一個示例中,sex和和age均被編碼為數字變數。如果存在字元變數,則必須使用model.matrix

格式化來自crr的結果

或當前crr不支援的輸出。 

黑色素瘤資料中的競爭風險迴歸-因果分析

刪失所有沒有引起關注的物件,在這種情況下是由於黑色素瘤死亡,並且照常使用coxph。因此,現在對因其他原因死亡的患者進行鍼對特定原因的風險評估方法以應對競爭風險。

第4部分:高階主題

 涵蓋的內容

  • 生存分析的基礎知識,包括Kaplan-Meier生存函式和Cox迴歸

  • 地標分析和時間相關協變數

  • 競爭風險分析的累積發生率和迴歸

還有什麼?

可能會出現很多零碎的東西 :

  1. 評估比例風險假設

  2. 生存率繪製平滑的生存圖XX

  3. 有條件的生存

評估比例風險

Cox比例風險迴歸模型的一個假設是,在整個隨訪過程中,風險在每個時間點都是成比例的。我們如何檢查資料是否符合此假設?

使用cox.zph生存包中的功能。結果有兩點:

  1. 每個協變數的效果是否隨時間變化的假設檢驗,以及一次所有協變數的全域性檢驗。

  • 這是通過證明協變數和log(time)之間的互動作用來完成的

  • 顯著的p值表示違反了比例風險假設

  1. Schoenfeld殘差圖

  • 偏離零坡度線的證據表明違反了比例風險假設

print(cz)
  1.   ##            rho chisq     p
  2.   ## sex     0.1236 2.452 0.117
  3.   ## age    -0.0275 0.129 0.719
  4.   ## GLOBAL      NA 2.651 0.266
  5.   ``````
  6.   plot(cz)

平滑的生存圖-生存分位數

有時可能想根據連續變數來視覺化生存估計。求 生存資料的分位數。預設分位數是p = 0.5中位生存期。

  • x代表事件

  • o代表刪失

  • 該線是根據年齡的平均存活率的平滑估計

條件生存

有時,在已經存活了一段時間的患者中產生存活率估計值很有意義。

Zabor, E., Gonen, M., Chapman, P., & Panageas, K. (2013). Dynamic prognostication using conditional survival estimates. Cancer, 119(20), 3589-3592.

條件生存估計

讓我們將生存期定為6個月

  1.   map_df(
  2.     prob_times, 
  3.     ~conditional\_surv\_est(
  4.       basekm = fit1, 
  5.       t1 = 182.625, 
  6.       t2 = .x) 
  7.     ) %>% 
  8.     mutate(months = round(prob_times / 30.4)) %>% 
  9.     select(months, everything()) %>% 
  10.     kable()

條件生存圖

我們還可以根據不同的生存時間長度視覺化條件生存資料。 

所得出的曲線在我們每次進行條件調整時都有一條生存曲線。在這種情況下,第一條線是總體生存曲線,因為它是根據時間0進行調節的。

資料獲取

在下面公眾號後臺回覆“肺癌資料”,可獲取完整資料。