【視頻】R語言生存分析原理與晚期肺癌患者分析案例|數(shù)據(jù)分享|附代碼數(shù)據(jù)
原文鏈接:http://tecdat.cn/?p=10278
最近我們被客戶要求撰寫關于生存分析的研究報告,包括一些圖形和統(tǒng)計輸出。
生存分析(也稱為工程中的可靠性分析)的目標是在協(xié)變量和事件時間之間建立聯(lián)系
生存分析的名稱源于臨床研究,其中預測死亡時間,即生存,通常是主要目標。
生存分析是一種回歸問題(人們想要預測一個連續(xù)值),但有一個轉折點。它與傳統(tǒng)回歸的不同之處在于,在生存分析中,結果變量既有一個事件,也有一個與之相關的時間值,部分訓練數(shù)據(jù)只能被部分觀察——它們是被刪失的。本文用R語言生存分析晚期肺癌患者數(shù)據(jù)
普通最小二乘回歸方法不足,因為事件發(fā)生的時間通常不是正態(tài)分布的,并且模型無法處理刪失,但這在生存數(shù)據(jù)中很常見。
為什么要做生存分析:右刪失
在某些情況下,可能無法觀察到事件時間:這通常稱為?右刪失。在以死亡為事件的臨床試驗中,當發(fā)生以下情況之一時,就會發(fā)生這種情況。1。當一定數(shù)量的參與者死亡時,研究結束。2。參與者退出研究。3。?研究達到預定的結束時間,并且一些參與者存活到結束。在每種情況下,幸存的參與者離開研究后,我們都不知道他們會發(fā)生什么。然后我們有一個問題:
當對于某些個體,我們只觀察到他們的事件時間的下限時,我們?nèi)绾螌?jīng)驗分布進行建?;蜻M行非負回歸?
上圖說明了右刪失。對于參與者 1,我們看到他們何時死亡。參與者 2 退出了,我們知道他們一直活到那時,但不知道后來發(fā)生了什么。對于參與者 3,我們知道他們活到了預定的研究結束,但又不知道之后發(fā)生了什么。
生存函數(shù)和風險函數(shù)
生存分析中的兩個關鍵工具是生存函數(shù)和風險函數(shù)。
生存函數(shù):它是一個函數(shù),用于給出我們有興趣知道的任何對象是否會在任何指定時間之后存活的概率。在數(shù)學上它可以由以下公式表示?
其中 S(t) 是一個生存函數(shù),其中 T 是一個連續(xù)隨機變量,是一個事件的時間。F(t) 是區(qū)間[0,∞) 上的累積分布函數(shù)。
我們也可以用風險函數(shù)來寫生存函數(shù)。假設事件尚未發(fā)生?,風險率λ(t) 是事件在時間t發(fā)生的瞬時概率的主要值。
那么關鍵問題是如何估計風險和/或生存函數(shù)。
Kaplan Meier的非參數(shù)估計
在非參數(shù)生存分析中,我們要估計生存函數(shù)沒有協(xié)變量,并且有刪失。如果我們沒有刪失,我們可以從經(jīng)驗 CDF 開始. 這個等式簡潔地表示:
有多少人隨著時間的推移而死亡? 那么生存函數(shù)就是:還有多少人還活著?
但是,我們無法回答一些人被時間t刪失時提出的這個問題.
雖然我們不一定知道有多少人在任意時間t幸存下來,我們知道研究中有多少人仍然處于風險之中。我們可以使用它來代替。將學習時間劃分區(qū)間, 其中每個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
library(survival)
什么是生存數(shù)據(jù)?
事件時間數(shù)據(jù)由不同的開始時間和結束時間組成。
癌癥的例子
從手術到死亡的時間
從治療開始到進展的時間
從響應到復發(fā)的時間
其他領域的例子
事件發(fā)生時間數(shù)據(jù)在許多領域都很常見,包括但不限于
從艾滋病毒感染到艾滋病發(fā)展的時間
心臟病發(fā)作的時間
藥物濫用發(fā)生的時間
機器故障時間
生存分析別名
由于生存分析在許多其他領域很常見,因此也有其他名稱
可靠性分析
持續(xù)時間分析
事件歷史分析
事件發(fā)生時間分析
肺數(shù)據(jù)集
數(shù)據(jù)包含來自北中部癌癥治療組的晚期肺癌患者。今天我們將用來演示方法的一些變量包括
時間:以天為單位的生存時間
狀態(tài):刪失狀態(tài)1 =刪失,2 =失效
性別:男= 1女= 2
刪失類型
某個主題可能由于以下原因而被刪失:
后續(xù)損失
退出研究
固定學習期結束前沒有活動
具體來說,這些是刪失的示例。
分配隨訪時間
受刪失的主題仍會提供信息,因此必須適當?shù)匕诜治鲋?/p>
隨訪時間的分布存在偏差,在接受檢查的患者和有事件的患者之間可能有所不同
生存數(shù)據(jù)的組成部分
對于主題ii:
活動時間Ti
刪失時間Ci
事件指標δi:
1,如果觀察到的事件(即??Ti≤CiTi≤Ci)
如果檢查,則為0(即??Ti>CiTi>Ci)
觀測時間Yi=min(Ti,Ci)Yi=min(Ti,Ci)
lung
數(shù)據(jù)中提供了觀察時間和事件指示
時間:以天為單位的生存時間(YiYi)
狀態(tài):刪失狀態(tài)1 =刪失,2 =死亡(δiδi)
在R中處理日期
數(shù)據(jù)通常帶有開始日期和結束日期,而不是預先計算的生存時間。第一步是確保將這些格式設置為R中的日期。
讓我們創(chuàng)建一個小的示例數(shù)據(jù)集,其中sx_date
包含手術日期和last_fup_date
上次隨訪日期的變量。
date_ex?<-?
??tibble(
????sx_date?=?c("2007-06-22",?"2004-02-13",?"2010-10-27"),?
????last_fup_date?=?c("2017-04-15",?"2018-07-04",?"2016-10-31")
????)
date_ex
##?#?A?tibble:?3?x?2##???sx_date????last_fup_date##???<chr>??????<chr>????????##?1?2007-06-22?2017-04-15???##?2?2004-02-13?2018-07-04???##?3?2010-10-27?2016-10-31
我們看到它們都是字符變量,通常都是這種情況,但是我們需要將它們格式化為日期。
格式化日期-基數(shù)R
date_ex?%>%?
??mutate(
????sx_date?=?as.Date(sx_date,?format?=?"%Y-%m-%d"),?
????last_fup_date?=?as.Date(last_fup_date,?format?=?"%Y-%m-%d")?
????)
##?#?A?tibble:?3?x?2##???sx_date????last_fup_date##???<date>?????<date>???????##?1?2007-06-22?2017-04-15???##?2?2004-02-13?2018-07-04???##?3?2010-10-27?2016-10-31
請注意,
R
格式必須包含分隔符和符號。例如,如果您的日期格式為m / d / Y,則需要format = "%m/%d/%Y"
格式化日期-lubridate程序包
我們還可以使用該lubridate
包來格式化日期。在這種情況下,請使用ymd
功能
date_ex?%>%?
??mutate(
????sx_date?=?ymd(sx_date),?
????last_fup_date?=?ymd(last_fup_date)
????)
##?#?A?tibble:?3?x?2##???sx_date????last_fup_date##???<date>?????<date>???????##?1?2007-06-22?2017-04-15???##?2?2004-02-13?2018-07-04???##?3?2010-10-27?2016-10-31
請注意,與基本
R
選項不同,不需要指定分隔符
計算生存時間?
現(xiàn)在日期已格式化,我們需要以某些單位(通常是幾個月或幾年)計算開始時間和結束時間之間的差。在base中R
,用于difftime
計算兩個日期之間的天數(shù),然后使用將其轉換為數(shù)字值as.numeric
。然后將除以365.25
年的平均天數(shù)轉換為年。
date_ex?%>%?
??mutate(
????os_yrs?=?
??????as.numeric(
????????difftime(last_fup_date,?
?????????????????sx_date,?
?????????????????units?=?"days"))?/?365.25????)
##?#?A?tibble:?3?x?3##???sx_date????last_fup_date?os_yrs##???<date>?????<date>?????????<dbl>##?1?2007-06-22?2017-04-15??????9.82##?2?2004-02-13?2018-07-04?????14.4?##?3?2010-10-27?2016-10-31??????6.01
計算生存時間?
操作員可以%--%
指定一個時間間隔,然后使用將該時間間隔轉換為經(jīng)過的秒數(shù)as.duration
,最后除以dyears(1)
,將其轉換為年數(shù),從而得出一年中的秒數(shù)。
##?#?A?tibble:?3?x?3##???sx_date????last_fup_date?os_yrs##???<date>?????<date>?????????<dbl>##?1?2007-06-22?2017-04-15??????9.82##?2?2004-02-13?2018-07-04?????14.4?##?3?2010-10-27?2016-10-31??????6.02
事件標標
對于生存數(shù)據(jù)的組成部分,我提到了事件指示器:
事件指標δiδi:
1,如果觀察到的事件(即??Ti≤CiTi≤Ci)
如果檢查,則為0(即??Ti>CiTi>Ci)
在lung
數(shù)據(jù)中,我們有:
狀態(tài):刪失狀態(tài)1 =刪失,2 =失效
生存函數(shù)
受試者可以存活超過指定時間的概率
S(t)=Pr(T>t)=1?F(t)S(t)=Pr(T>t)=1?F(t)
S(t)S(t):生存函數(shù)F(t)=Pr(T≤t)F(t)=Pr(T≤t):累積分布函數(shù)
理論上,生存函數(shù)是平滑的;在實踐中,我們以離散的時間尺度觀察事件。
生存概率
生存概率在某個時間,S(t)S(t),是存活超過該時間,考慮到個體已存活剛剛在此之前,時間的條件概率。
可以估計為當時活著但沒有損失的隨訪患者人數(shù)除以當時的活著患者人數(shù)
生存概率的Kaplan-Meier估計是這些條件概率的乘積
在時間0,生存概率為1,即??S(t0)=1S(t0)=1
創(chuàng)建生存對象
Kaplan-Meier方法是估計生存時間和概率的最常用方法。這是一種非參數(shù)方法,可產(chǎn)生階躍函數(shù),每次事件發(fā)生時,階躍下降。
創(chuàng)建一個生存對象。對于每個主題,將有一個條目作為生存時間,
+
如果主題是經(jīng)過刪失的,則后面跟一個。讓我們看一下前10個觀察值:
##??[1]??306???455??1010+??210???883??1022+??310???361???218???166
用Kaplan-Meier方法估算生存曲線
該
survfit
函數(shù)根據(jù)公式創(chuàng)建生存曲線。讓我們?yōu)檎麄€同類群組生成總體生存曲線,將其分配給object?f1
,然后查看names
該對象的:
names(f1)
##??[1]?"n"??????????"time"???????"n.risk"?????"n.event"????"n.censor"??
##??[6]?"surv"???????"std.err"????"cumhaz"?????"std.chaz"???"start.time"##?[11]?"type"???????"logse"??????"conf.int"???"conf.type"??"lower"?????
##?[16]?"upper"??????"call"
該survfit
對象將用于創(chuàng)建生存曲線的一些關鍵組件包括:
time
,其中包含每個時間間隔的起點和終點surv
,其中包含每個對應的生存概率?time
Kaplan-Meier圖?
現(xiàn)在, 繪制對象 獲得Kaplan-Meier圖。
plot(survfit(Surv(time,?status)?~?1,?data?=?lung),?
????
基數(shù)
R
中的默認圖顯示了具有相關置信區(qū)間(虛線)的階躍函數(shù)(實線)水平線代表間隔的生存時間
時間間隔由事件終止
垂直線的高度顯示累積概率的變化
帶有刻度線的經(jīng)過刪失的觀察結果會減少間隔之間的累積生存期。
Kaplan-Meier圖?
建立在上ggplot2
,并可用于創(chuàng)建Kaplan-Meier圖。
點擊標題查閱往期內(nèi)容
R語言生存分析數(shù)據(jù)分析可視化案例
左右滑動查看更多
01
02
03
04
默認圖??帶相關置信帶(陰影區(qū)域)的階躍函數(shù)(實線)。
默認情況下,顯示了被檢查患者的刻度線,在此示例中,該刻度線本身有些模糊,可以使用選項將其取消?
censor = FALSE
估計xx年生存
生存分析中經(jīng)常需要關注的一個數(shù)量是生存超過一定數(shù)量(xx)年的概率。
例如,要估算生存到11年的可能性
##?Call:?survfit(formula?=?Surv(time,?status)?~?1,?data?=?lung)##?##??time?n.risk?n.event?survival?std.err?lower?95%?CI?upper?95%?CI##???365?????65?????121????0.409??0.0358????????0.345????????0.486
我們發(fā)現(xiàn)本研究中11年生存的機率是41%。
同時顯示95%置信區(qū)間的相關上下限。
xx年生存率和生存曲線
11年存活率概率為在y軸上的點對應于11一年x軸的生存曲線。
Xx年生存率常常被錯誤估計
如果 使用“天真”的估計會怎樣?
228名患者中的121名到1年時死亡,因此:
-當 忽略42名患者在1年之前受到檢查的事實時, 會錯誤估計1個1個年生存率。
正確的估計生存概率-年為41%。
忽略刪失對xx年生存的影響
想象兩個研究,每個研究228個主題。每個研究中有165人死亡。一個沒有檢查(紅色線),63個病人被另一個(藍色線)檢查
忽略刪失會導致總體生存概率被高估,因為被刪失的受試者僅在部分隨訪時間內(nèi)提供信息,然后落入風險范圍之外,從而降低了生存的累積概率
估計中位生存時間
生存分析中經(jīng)常需要關注的另一個數(shù)量是平均生存時間,我們使用中位數(shù)對其進行量化。預計生存時間不會呈正態(tài)分布,因此平均值不是適當?shù)目偨Y。
##?Call:?survfit(formula?=?Surv(time,?status)?~?1,?data?=?lung)##?##???????n??events??median?0.95LCL?0.95UCL?##?????228?????165?????310?????285?????363
我們看到中位生存時間為310天。還會顯示95%置信區(qū)間的上限和下限。
中位生存時間和生存曲線
中位生存時間是生存概率為0.50?
中位生存率常常被錯誤估計
總結165例死亡患者的中位生存時間
##???median_surv##?1?????????226
當 忽略被檢查患者也有助于隨訪的事實時, 會得出錯誤的估計中值生存時間226天。
正確的中位生存時間估計是310天。
忽略刪失對中位數(shù)生存率的影響
忽略刪失會造成人為降低的生存曲線,因為排除了受刪失患者貢獻的隨訪時間(紫色線)
數(shù)據(jù)的真實生存曲線以
lung
藍色顯示,以進行比較
比較各組之間的生存時間
我們可以使用對數(shù)秩檢驗進行組間重要性檢驗
對數(shù)秩檢驗在整個隨訪時間內(nèi)平均權衡觀察結果,是比較組間生存時間的最常用方法
根據(jù)研究問題,有些版本可能會更重視早期或后期的隨訪,可能更合適
我們使用 函數(shù)獲得對數(shù)秩p值。例如,我們可以根據(jù)lung
數(shù)據(jù)中的性別測試是否存在生存時間差異
##?Call:##?##?????????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.001
從survdiff對象中提取信息
從 結果中提取p值
1?-?pchisq(sd$chisq,?length(sd$n)?-?1)
##?[1]?0.001311165
返回格式化的p值
##?[1]?0.001
Cox回歸模型
我們可能想量化單個變量的效應大小,或者將多個變量包括在回歸模型中以說明多個變量的效應。
Cox回歸模型是半?yún)?shù)模型,可用于擬合具有生存結果的單變量和多變量回歸模型。
h(t)h(t):危險或事件發(fā)生的瞬時速率h0(t)h0(t):基本基準危險
該模型的一些關鍵假設:
非信息刪失
比例危險
注意:也可以使用用于生存結果的參數(shù)回歸模型,但是本培訓將不涉及這些模型。
我們可以使用coxph
函數(shù)擬合生存數(shù)據(jù)的回歸模型,該函數(shù)Surv
在左側使用一個對象,而在右側具有用于回歸公式的標準語法R
。
##?Call:##?##????????coef?exp(coef)?se(coef)??????z???????p##?sex?-0.5310????0.5880???0.1672?-3.176?0.00149##?##?Likelihood?ratio?test=10.63??on?1?df,?p=0.001111##?n=?228,?number?of?events=?165
格式化Cox回歸結果
我們可以看到輸出的整潔版本broom
:
或使用
危險比
來自Cox回歸模型的關注數(shù)量是危險比(HR)。HR表示在任何特定時間點兩組之間的危險比率。
HR被解釋為感興趣事件中那些仍處于事件風險中的事件的瞬時發(fā)生率。
如果您有一個回歸參數(shù)ββ(來自
estimate
我們的列coxph
),則HR =?經(jīng)驗值(β)經(jīng)驗值?(β)。HR <1表示死亡危險降低,而HR> 1表示死亡危險增加。
因此,我們的HR = 0.59意味著在任何給定時間,女性死亡的人數(shù)大約是男性的0.6倍。
第2部分:地標分析和時間相關協(xié)變量
在第1部分中,我們介紹了使用對數(shù)秩檢驗和Cox回歸來檢驗感興趣的協(xié)變量與生存結果之間的關聯(lián)。
示例:腫瘤反應
示例:從治療開始就測量總生存期,關注的是對治療的完全反應與生存之間的關聯(lián)。
Anderson等人(JCO,1983)描述了在這種情況下,為什么傳統(tǒng)方法(如對數(shù)秩檢驗或Cox回歸)偏向于響應者,并提出了劃時代的方法。
界標方法中的零假設是,從界標生存的過程不依賴于界標的響應狀態(tài)。
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.
其他例子
癌癥研究中可能尚未關注的其他一些可能的協(xié)變量包括:
移植失敗
移植物抗宿主病
第二次切除
輔助治療
合規(guī)
不良事件
示例數(shù)據(jù)?
137例骨髓移植患者的數(shù)據(jù)。變量包括:
T1
?死亡時間或最后一次隨訪時間(天)delta1
死亡指標;1死0活TA
?急性移植物抗宿主病的時間(以天為單位)deltaA
急性移植物抗宿主病指標;1-發(fā)展為急性移植物抗宿主病,0-從未發(fā)展為急性移植物抗宿主病
讓我們加載數(shù)據(jù)以供整個示例使用
地標法
選擇基線之后的固定時間作為界標時間。注意:應在檢查數(shù)據(jù)之前根據(jù)臨床信息進行操作
那些人群的子集至少跟蹤到里程碑時間。注意:請務必在地標時間之前報告由于關注或刪失事件而排除的號碼。
計算具有里程碑意義的時間,并應用傳統(tǒng)的對數(shù)秩檢驗或Cox回歸
在BMT
數(shù)據(jù)感興趣的是急性移植物抗宿主?。╝GVHD)和存活之間的關聯(lián)。但是aGVHD是在移植后進行評估的,這是我們的基線,也就是后續(xù)隨訪的開始時間。
步驟1選擇地標時間
通常,aGVHD發(fā)生在移植后的前90天內(nèi),因此我們使用90天的界標。
人們對急性移植物抗宿主?。╝GVHD)與生存之間的關系感興趣。但是aGVHD是在移植后進行評估的,這是我們的基線,也就是后續(xù)隨訪的開始時間。
第2步:至少跟蹤到里程碑時間之前的人群的子集
這將我們的樣本量從137減少到122。
所有15位被排除的患者均在90天里程碑之前死亡
人們對急性移植物抗宿主?。╝GVHD)與生存之間的關系感興趣。但是aGVHD是在移植后進行評估的,這是我們的基線,也就是后續(xù)隨訪的開始時間。
步驟3根據(jù)地標計算隨訪時間,并應用傳統(tǒng)方法。
使用BMT數(shù)據(jù)的Cox回歸界標示例
在Cox回歸中, 可以使用中的subset
選項coxph
來排除那些在標志性時間內(nèi)沒有被隨訪的患者
時間相關協(xié)變量
界標分析的替代方法是合并時間相關的協(xié)變量。這可能更適合
協(xié)變量的值隨時間變化
沒有明顯的里程碑時間
時間相關協(xié)變量數(shù)據(jù)設置
對時間相關協(xié)變量的分析R
需要建立特殊的數(shù)據(jù)集。?
BMT
數(shù)據(jù)中沒有ID變量,這是創(chuàng)建特殊數(shù)據(jù)集所必需的,因此請創(chuàng)建一個名為的變量my_id
。
將tmerge
函數(shù)與event
和函數(shù)一起使用tdc
可創(chuàng)建特殊數(shù)據(jù)集。
tmerge
?為每個患者的不同協(xié)變量值創(chuàng)建一個具有多個時間間隔的長數(shù)據(jù)集event
?創(chuàng)建新的事件指示器,以與新創(chuàng)建的時間間隔一致tdc
?創(chuàng)建與時間相關的協(xié)變量指標,以與新創(chuàng)建的時間間隔一致
時間相關協(xié)變量-單例患者
要了解其作用,讓我們看一下前5名患者的數(shù)據(jù)。
?##???my_id???T1?delta1???TA?deltaA##?1?????1?2081??????0???67??????1##?2?????2?1602??????0?1602??????0##?3?????3?1496??????0?1496??????0##?4?????4?1462??????0???70??????1##?5?????5?1433??????0?1433??????0
這些相同患者的新數(shù)據(jù)集
##???my_id???T1?delta1?id?tstart?tstop?death?agvhd##?1?????1?2081??????0??1??????0????67?????0?????0##?2?????1?2081??????0??1?????67??2081?????0?????1##?3?????2?1602??????0??2??????0??1602?????0?????0##?4?????3?1496??????0??3??????0??1496?????0?????0##?5?????4?1462??????0??4??????0????70?????0?????0##?6?????4?1462??????0??4?????70??1462?????0?????1##?7?????5?1433??????0??5??????0??1433?????0?????0
時間相關協(xié)變量-Cox回歸
現(xiàn)在,我們可以分析這個時間依賴性協(xié)照常使用Cox回歸與coxph?
摘要
我們發(fā)現(xiàn),使用標志性分析或時間依賴性協(xié)變量,急性移植物抗宿主病與死亡無顯著相關性。
通常,人們會希望使用地標分析對單個協(xié)變量進行可視化, 使用帶有時間相關協(xié)變量的Cox回歸進行單變量和多變量建模。
第3部分:競爭風險
什么是競爭風險?
當對象在事件發(fā)生時間設置中發(fā)生多個可能的事件時
例子:
復發(fā)
因疾病死亡
因其他原因死亡
治療反應
在任何給定的研究中,所有這些(或其中一些 以及其他)可能都是可能的事件。
所以有什么問題?
事件時間之間未觀察到的依賴性是導致需要特殊考慮的基本問題。
例如,可以想象復發(fā)的患者更有可能死亡,因此復發(fā)時間和死亡時間將不是獨立事件。
競爭風險的背景
存在多種潛在結果時的兩種分析方法:
給定事件的特定于原因的危險:這表示未因其他事件而失敗的事件中事件的每單位時間的發(fā)生率
給定事件的累積發(fā)生率:這表示事件每單位時間的發(fā)生率以及競爭事件的影響
這些方法中的每一種都可能僅闡明數(shù)據(jù)的一個重要方面,而有可能使其他方面難以理解,因此所選的方法應取決于感興趣的問題。
?黑色素瘤數(shù)據(jù)示例
它包含變量:
time
?生存時間以天為單位,可能經(jīng)過刪失。status
?1例死于黑色素瘤,2例存活,3例因其他原因死亡。sex
?1 =男性,0 =女性。age
?年歲。year
?操作。thickness
?腫瘤厚度(毫米)。ulcer
?1 =存在,0 =不存在。
黑色素瘤數(shù)據(jù)的累積發(fā)生率
在競爭風險的背景下估算累積發(fā)生率。
##?Estimates?and?Variances:##?$est##???????????1000???????2000???????3000??????4000??????5000##?1?1?0.12745714?0.23013963?0.30962017?0.3387175?0.3387175##?1?3?0.03426709?0.05045644?0.05811143?0.1059471?0.1059471##?##?$var##?????????????1000?????????2000?????????3000????????4000????????5000##?1?1?0.0005481186?0.0009001172?0.0013789328?0.001690760?0.001690760##?1?3?0.0001628354?0.0002451319?0.0002998642?0.001040155?0.001040155
繪制累積發(fā)生率-基數(shù)R
生成 默認值的基本圖。
plot(ci_fit)
繪制累積發(fā)生率?
比較組之間的累積發(fā)生率
用于組間測試。
例如,Melanoma
根據(jù)ulcer
潰瘍的存在與否比較結果。測試結果可以在中找到Tests
。
ci_ulcer[["Tests"]]
##????????stat???????????pv?df##?1?26.120719?3.207240e-07??1##?3??0.158662?6.903913e-01??1
按組繪制累積發(fā)生率?
按組繪制累積發(fā)生率-手動
_請注意,_我個人發(fā)現(xiàn)該ggcompetingrisks
功能缺少自定義功能,尤其是與相比ggsurvplot
。我通常會自己做圖,首先創(chuàng)建cuminc
擬合結果的整潔數(shù)據(jù)集,然后再繪制結果。有關底層代碼的詳細信息,請參見此演示文稿的
繪制單個事件類型?
通常,只有一種類型的事件會引起人們的興趣,盡管我們?nèi)砸紤]競爭事件。在那種情況下,感興趣的事件可以單獨繪制。同樣,我首先通過創(chuàng)建cuminc
擬合結果的整潔數(shù)據(jù)集,然后繪制結果來手動執(zhí)行此操作。有關底層代碼的詳細信息,請參見此演示文稿的源代碼。
在風險表中添加數(shù)字
您可能想將風險表的數(shù)量添加到累積發(fā)生率圖中,而據(jù)我所知,沒有簡單的方法可以做到這一點。請參閱此演示文稿的源代碼中的一個示例
競爭風險回歸
兩種方法:
特定原因風險
當前沒有事件的受試者中給定事件類型的瞬時發(fā)生率
使用Cox回歸估算
Subdistribution子分布風險
給定類型事件在沒有經(jīng)歷過此類事件的受試者中的瞬時發(fā)生率
使用Fine-Gray回歸估算
黑色素瘤數(shù)據(jù)中的競爭風險回歸-子分布風險法Subdistribution
假設我們有興趣研究年齡和性別對黑色素瘤死亡的影響,而其他原因的死亡則是競爭事件。
crr
?需要指定協(xié)變量作為矩陣如果 關注多個事件,則可以使用
failcode
選項請求其他事件的結果,默認情況下會返回failcode = 1
shr_fit
##?convergence:??TRUE?##?coefficients:##?????sex?????age?##?0.58840?0.01259?##?standard?errors:##?[1]?0.271800?0.009301##?two-sided?p-values:##??sex??age?##?0.03?0.18
在上一個示例中,sex
和和age
均被編碼為數(shù)字變量。如果存在字符變量,則必須使用model.matrix
格式化來自crr的結果
或當前crr
不支持的輸出。?
黑色素瘤數(shù)據(jù)中的競爭風險回歸-因果分析
刪失所有沒有引起關注的對象,在這種情況下是由于黑色素瘤死亡,并且照常使用coxph
。因此,現(xiàn)在對因其他原因死亡的患者進行針對特定原因的風險評估方法以應對競爭風險。
第4部分:高級主題
?涵蓋的內(nèi)容
生存分析的基礎知識,包括Kaplan-Meier生存函數(shù)和Cox回歸
地標分析和時間相關協(xié)變量
競爭風險分析的累積發(fā)生率和回歸
還有什么?
可能會出現(xiàn)很多零碎的東西 :
評估比例風險假設
生存率繪制平滑的生存圖XX
有條件的生存
評估比例風險
Cox比例風險回歸模型的一個假設是,在整個隨訪過程中,風險在每個時間點都是成比例的。我們?nèi)绾螜z查數(shù)據(jù)是否符合此假設?
使用cox.zph
生存包中的功能。結果有兩點:
每個協(xié)變量的效果是否隨時間變化的假設檢驗,以及一次所有協(xié)變量的全局檢驗。
這是通過證明協(xié)變量和log(time)之間的交互作用來完成的
顯著的p值表示違反了比例風險假設
Schoenfeld殘差圖
偏離零坡度線的證據(jù)表明違反了比例風險假設
print(cz)
##????????????rho?chisq?????p##?sex?????0.1236?2.452?0.117##?age????-0.0275?0.129?0.719##?GLOBAL??????NA?2.651?0.266``````
plot(cz)
平滑的生存圖-生存分位數(shù)
有時可能想根據(jù)連續(xù)變量來可視化生存估計。求 生存數(shù)據(jù)的分位數(shù)。默認分位數(shù)是p = 0.5
中位生存期。
x代表事件
o代表刪失
該線是根據(jù)年齡的平均存活率的平滑估計
條件生存
有時,在已經(jīng)存活了一段時間的患者中產(chǎn)生存活率估計值很有意義。
Zabor, E., Gonen, M., Chapman, P., & Panageas, K. (2013). Dynamic prognostication using conditional survival estimates. Cancer, 119(20), 3589-3592.
條件生存估計
讓我們將生存期定為6個月
map_df(
??prob_times,?
??~conditional_surv_est(
????basekm?=?fit1,?
????t1?=?182.625,?
????t2?=?.x)?
??)?%>%?
??mutate(months?=?round(prob_times?/?30.4))?%>%?
??select(months,?everything())?%>%?
??kable()
條件生存圖
我們還可以根據(jù)不同的生存時間長度可視化條件生存數(shù)據(jù)。?
所得出的曲線在我們每次進行條件調(diào)整時都有一條生存曲線。在這種情況下,第一條線是總體生存曲線,因為它是根據(jù)時間0進行調(diào)節(jié)的。
數(shù)據(jù)獲取
在下面公眾號后臺回復“肺癌****數(shù)據(jù)”,可獲取完整數(shù)據(jù)。
本文摘選?《?R語言中的生存分析Survival analysis晚期肺癌患者?》?,點擊“閱讀原文”獲取全文完整資料。
點擊標題查閱往期內(nèi)容
R語言使用限制平均生存時間RMST比較兩條生存曲線分析肝硬化患者
生存分析模型的時間依賴性ROC曲線可視化R語言生存分析: 時變競爭風險模型分析淋巴瘤患者
R語言生存分析可視化分析
R語言中生存分析模型的時間依賴性ROC曲線可視化
R語言生存分析數(shù)據(jù)分析可視化案例
R語言ggsurvplot繪制生存曲線報錯 : object of type ‘symbol‘ is not subsettab
R語言如何在生存分析與Cox回歸中計算IDI,NRI指標
R語言繪制生存曲線估計|生存分析|如何R作生存曲線圖
R語言解釋生存分析中危險率和風險率的變化
R語言中的生存分析Survival analysis晚期肺癌患者4例