原文連結:http://tecdat.cn/?p=10278
最近我們被客戶要求撰寫關於生存分析的研究報告,包括一些圖形和統計輸出。
生存分析(也稱為工程中的可靠性分析)的目標是在協變量和事件時間之間建立聯繫
生存分析的名稱源於臨床研究,其中預測死亡時間,即生存,通常是主要目標。
視頻:R語言生存分析原理與晚期肺癌患者分析案例
**,時長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 估計示例。
生存分析用於各種領域
例如:
在癌症研究中,典型的研究問題如下:
本演示文稿將介紹生存分析 ,參考:
Clark, T., Bradburn, M., Love, S., & Altman, D. (2003). Survival analysis part I: Basic concepts and first analyses. 232-238. ISSN 0007-0920.
我們今天將使用的一些軟體包包括:
library(survival)
什麼是生存數據?
事件時間數據由不同的開始時間和結束時間組成。
癌症的例子
其他領域的例子
事件發生時間數據在許多領域都很常見,包括但不限於
生存分析別名
由於生存分析在許多其他領域很常見,因此也有其他名稱
肺數據集
數據包含來自北中部癌症治療組的晚期肺癌患者。今天我們將用來演示方法的一些變量包括
刪失類型
某個主題可能由於以下原因而被刪失:
具體來說,這些是刪失的示例。
分配隨訪時間
生存數據的組成部分
對於主題ii:
lung數據中提供了觀察時間和事件指示
在R中處理日期
數據通常帶有開始日期和結束日期,而不是預先計算的生存時間。第一步是確保將這些格式設置為R中的日期。
讓我們創建一個小的示例數據集,其中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
##
## 1 2007-06-22 2017-04-15
## 2 2004-02-13 2018-07-04
## 3 2010-10-27 2016-10-31
我們看到它們都是字符變量,通常都是這種情況,但是我們需要將它們格式化為日期。
格式化日期-基數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
##
## 1 2007-06-22 2017-04-15
## 2 2004-02-13 2018-07-04
## 3 2010-10-27 2016-10-31
格式化日期-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
##
## 1 2007-06-22 2017-04-15
## 2 2004-02-13 2018-07-04
## 3 2010-10-27 2016-10-31
計算生存時間
現在日期已格式化,我們需要以某些單位(通常是幾個月或幾年)計算開始時間和結束時間之間的差。在base中R,用於difftime計算兩個日期之間的天數,然後使用將其轉換為數字值as.numeric。然後將除以365.25年的平均天數轉換為年。
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
##
## 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
計算生存時間
操作員可以%--%指定一個時間間隔,然後使用將該時間間隔轉換為經過的秒數as.duration,最後除以dyears(1),將其轉換為年數,從而得出一年中的秒數。
## # A tibble: 3 x 3
## sx_date last_fup_date os_yrs
##
## 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
事件標標
對於生存數據的組成部分,我提到了事件指示器:
事件指標δiδi:
在lung數據中,我們有:
生存函數
受試者可以存活超過指定時間的機率
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):累積分布函數
理論上,生存函數是平滑的;在實踐中,我們以離散的時間尺度觀察事件。
生存機率
創建生存對象
Kaplan-Meier方法是估計生存時間和機率的最常用方法。這是一種非參數方法,可產生階躍函數,每次事件發生時,階躍下降。
## [1] 306 455 1010+ 210 883 1022+ 310 361 218 166
用Kaplan-Meier方法估算生存曲線
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對象將用於創建生存曲線的一些關鍵組件包括:
Kaplan-Meier圖
現在, 繪製對象 獲得Kaplan-Meier圖。
plot(survfit(Surv(time, status) ~ 1, data = lung),
Kaplan-Meier圖
建立在上ggplot2,並可用於創建Kaplan-Meier圖。
點擊標題查閱往期內容
R語言生存分析數據分析可視化案例
左右滑動查看更多
01
02
03
04
估計xx年生存
生存分析中經常需要關注的一個數量是生存超過一定數量(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
我們發現本研究中11年生存的機率是41%。
同時顯示95%置信區間的相關上下限。
xx年生存率和生存曲線
11年存活率機率為在y軸上的點對應於11一年x軸的生存曲線。
Xx年生存率常常被錯誤估計
如果 使用「天真」的估計會怎樣?
228名患者中的121名到1年時死亡,因此:
-當 忽略42名患者在1年之前受到檢查的事實時, 會錯誤估計1個1個年生存率。
忽略刪失對xx年生存的影響
估計中位生存時間
生存分析中經常需要關注的另一個數量是平均生存時間,我們使用中位數對其進行量化。預計生存時間不會呈正態分布,因此平均值不是適當的總結。
## Call: survfit(formula = Surv(time, status) ~ 1, data = lung)
##
## n events median 0.95LCL 0.95UCL
## 228 165 310 285 363
我們看到中位生存時間為310天。還會顯示95%置信區間的上限和下限。
中位生存時間和生存曲線
中位生存時間是生存機率為0.50
中位生存率常常被錯誤估計
總結165例死亡患者的中位生存時間
## median_surv
## 1 226
忽略刪失對中位數生存率的影響
比較各組之間的生存時間
我們使用 函數獲得對數秩p值。例如,我們可以根據lung數據中的性別測試是否存在生存時間差異
## 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回歸模型是半參數模型,可用於擬合具有生存結果的單變量和多變量回歸模型。
h(t)h(t):危險或事件發生的瞬時速率h0(t)h0(t):基本基準危險
該模型的一些關鍵假設:
注意:也可以使用用於生存結果的參數回歸模型,但是本培訓將不涉及這些模型。
我們可以使用coxph函數擬合生存數據的回歸模型,該函數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:
或使用
危險比
在第1部分中,我們介紹了使用對數秩檢驗和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例骨髓移植患者的數據。變量包括:
讓我們加載數據以供整個示例使用
地標法
在BMT數據感興趣的是急性移植物抗宿主病(aGVHD)和存活之間的關聯。但是aGVHD是在移植後進行評估的,這是我們的基線,也就是後續隨訪的開始時間。
步驟1選擇地標時間
通常,aGVHD發生在移植後的前90天內,因此我們使用90天的界標。
人們對急性移植物抗宿主病(aGVHD)與生存之間的關係感興趣。但是aGVHD是在移植後進行評估的,這是我們的基線,也就是後續隨訪的開始時間。
第2步:至少跟蹤到里程碑時間之前的人群的子集
這將我們的樣本量從137減少到122。
人們對急性移植物抗宿主病(aGVHD)與生存之間的關係感興趣。但是aGVHD是在移植後進行評估的,這是我們的基線,也就是後續隨訪的開始時間。
步驟3根據地標計算隨訪時間,並應用傳統方法。
使用BMT數據的Cox回歸界標示例
在Cox回歸中, 可以使用中的subset選項coxph來排除那些在標誌性時間內沒有被隨訪的患者
時間相關協變量
界標分析的替代方法是合併時間相關的協變量。這可能更適合
時間相關協變量數據設置
對時間相關協變量的分析R需要建立特殊的數據集。
BMT數據中沒有ID變量,這是創建特殊數據集所必需的,因此請創建一個名為的變量my_id。
將tmerge函數與event和函數一起使用tdc可創建特殊數據集。
時間相關協變量-單例患者
要了解其作用,讓我們看一下前5名患者的數據。
## 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
這些相同患者的新數據集
## 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
時間相關協變量-Cox回歸
現在,我們可以分析這個時間依賴性協照常使用Cox回歸與coxph
摘要
我們發現,使用標誌性分析或時間依賴性協變量,急性移植物抗宿主病與死亡無顯著相關性。
通常,人們會希望使用地標分析對單個協變量進行可視化, 使用帶有時間相關協變量的Cox回歸進行單變量和多變量建模。
什麼是競爭風險?
當對象在事件發生時間設置中發生多個可能的事件時
例子:
在任何給定的研究中,所有這些(或其中一些 以及其他)可能都是可能的事件。
所以有什麼問題?
事件時間之間未觀察到的依賴性是導致需要特殊考慮的基本問題。
例如,可以想像復發的患者更有可能死亡,因此復發時間和死亡時間將不是獨立事件。
競爭風險的背景
存在多種潛在結果時的兩種分析方法:
這些方法中的每一種都可能僅闡明數據的一個重要方面,而有可能使其他方面難以理解,因此所選的方法應取決於感興趣的問題。
黑色素瘤數據示例
它包含變量:
黑色素瘤數據的累積發生率
在競爭風險的背景下估算累積發生率。
## 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
繪製累積發生率-基數R
生成 默認值的基本圖。
plot(ci_fit)
繪製累積發生率
比較組之間的累積發生率
用於組間測試。
例如,Melanoma根據ulcer潰瘍的存在與否比較結果。測試結果可以在中找到Tests。
ci_ulcer[["Tests"]]
## stat pv df
## 1 26.120719 3.207240e-07 1
## 3 0.158662 6.903913e-01 1
按組繪製累積發生率
按組繪製累積發生率-手動
_請注意,_我個人發現該ggcompetingrisks功能缺少自定義功能,尤其是與相比ggsurvplot。我通常會自己做圖,首先創建cuminc擬合結果的整潔數據集,然後再繪製結果。有關底層代碼的詳細信息,請參見此演示文稿的
繪製單個事件類型
通常,只有一種類型的事件會引起人們的興趣,儘管我們仍要考慮競爭事件。在那種情況下,感興趣的事件可以單獨繪製。同樣,我首先通過創建cuminc擬合結果的整潔數據集,然後繪製結果來手動執行此操作。有關底層代碼的詳細信息,請參見此演示文稿的原始碼。
在風險表中添加數字
您可能想將風險表的數量添加到累積發生率圖中,而據我所知,沒有簡單的方法可以做到這一點。請參閱此演示文稿的原始碼中的一個示例
競爭風險回歸
兩種方法:
黑色素瘤數據中的競爭風險回歸-子分布風險法Subdistribution
假設我們有興趣研究年齡和性別對黑色素瘤死亡的影響,而其他原因的死亡則是競爭事件。
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均被編碼為數字變量。如果存在字符變量,則必須使用model.matrix
格式化來自crr的結果
或當前crr不支持的輸出。
黑色素瘤數據中的競爭風險回歸-因果分析
刪失所有沒有引起關注的對象,在這種情況下是由於黑色素瘤死亡,並且照常使用coxph。因此,現在對因其他原因死亡的患者進行針對特定原因的風險評估方法以應對競爭風險。
涵蓋的內容
還有什麼?
可能會出現很多零碎的東西 :
評估比例風險
Cox比例風險回歸模型的一個假設是,在整個隨訪過程中,風險在每個時間點都是成比例的。我們如何檢查數據是否符合此假設?
使用cox.zph生存包中的功能。結果有兩點:
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)
平滑的生存圖-生存分位數
有時可能想根據連續變量來可視化生存估計。求 生存數據的分位數。默認分位數是p = 0.5中位生存期。
條件生存
有時,在已經存活了一段時間的患者中產生存活率估計值很有意義。
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()
條件生存圖
我們還可以根據不同的生存時間長度可視化條件生存數據。
所得出的曲線在我們每次進行條件調整時都有一條生存曲線。在這種情況下,第一條線是總體生存曲線,因為它是根據時間0進行調節的。
點擊標題查閱往期內容
R語言使用限制平均生存時間RMST比較兩條生存曲線分析肝硬化患者
生存分析模型的時間依賴性ROC曲線可視化R語言生存分析: 時變競爭風險模型分析淋巴瘤患者
R語言生存分析可視化分析
R語言中生存分析模型的時間依賴性ROC曲線可視化
R語言生存分析數據分析可視化案例
R語言ggsurvplot繪製生存曲線報錯 : object of type 『symbol『 is not subsettab
R語言如何在生存分析與Cox回歸中計算IDI,NRI指標
R語言繪製生存曲線估計|生存分析|如何R作生存曲線圖
R語言解釋生存分析中危險率和風險率的變化
R語言中的生存分析Survival analysis晚期肺癌患者4例