R語言做生存分析:繪製 Kaplan-Meirer生存曲線和實現Log-rank檢驗

2021-01-09 言爸說育兒

劉老師總結的統計教程合集，可以節省你不少時間

一般生存分析文章的統計部分通常會這麼寫：採用Kaplan-Meier法計算生存率和中位生存期，採用log-rank檢驗生存率差異，Cox回歸法進行多因素分析。檢驗水準α= 0.05。

我們首先看看什麼是生存分析？

因為無法在短時間內評價慢性病患者的預後，所以通常情況下不會簡單地採用治癒率、病死率等指標，而是對患者進行隨訪，分析一定的時間之後患者生存或死亡的情況，這種將事件的結果和出現這一結果所經歷的時間結合起來分析的方法，稱為生存分析 (Survival Analysis)。

其實，生存分析應用非常廣泛，可以分析在不同試驗條件下，研究對象「生存時間」的分布情況，進而了解試驗條件對生存時間的影響。這裡的「生存時間」不專指人或動物的生命持續時間，而是泛指某個事件發生前所經歷的時間。如一個工人從下崗後到實現再就業的時間；一臺汽車從開始使用到發生第一次故障的時間；一個病人從確診患病到死亡的時間等。

Kaplan-Meier法簡稱K-M法,又稱乘積極限法(Product-limit Estimate)是生存分析方法中最常用的一種,主要用於估計患者生存率和繪製生存曲線。Kaplan-Meier曲線,以生存時間為橫軸,生存率S (tk)為縱軸,繪製而成的連續型的階梯形曲線，用以說明生存時間與生存率之間的關係。生存曲線,一般是平滑而水平延伸的,當某個時間點一旦有患者發生終點事件(如死亡)，曲線就會垂直下降,下降幅度是該時間點上患者發生終點事件例數和上一個時間節點後隨訪的患者樣本量的比。

當兩組或多組生存曲線進行比較時,常用的假設檢驗方法是對數秩檢驗(log-rank test ) :又稱時序檢驗,屬於非參數檢驗,用於比較兩組或多組生存曲線或生存時間是否相同，檢驗統計量為卡方。

下面我們一起來看一下如何用R語言繪製Kaplan-Meier曲線和實現log-rank檢驗。

一、 安裝和加載R包

繪製Kaplan-Meier生存曲線需要安裝R包：survminer和survival。

install.packages("survminer")

install.packages("survival")

library(survminer)

library(survival)

二、導入內置數據集

我們使用survival包的lung數據集進行演示。

data(lung) # 加載lung數據集

View(lung) # 查看數據集

str(lung)

三、擬合生存曲線

3.1 創建生存對象

在survival包中先使用Surv()函數創建生存對象，生存對象是將事件時間和結局信息合併在一起的數據結構。

attach(lung) # 綁定數據集

Surv(time,status) # 創建生存對象

在上面輸出的生存對象中，帶"+"號的表示右刪失數據。

3.2 擬合曲線

R中使用survfit()函數來擬合生存曲線。

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

fit # 查看擬合曲線信息

還可以使用summary()函數輸出更多詳細信息。

summary(fit)

四、繪製Kaplan-Meier生存曲線

ggsurvplot(fit, data = lung)

ggsurvplot(fit, data = lung,

surv.median.line = "hv") # 增加中位生存時間

ggsurvplot(fit, data = lung,

surv.median.line = "hv", # 增加中位生存時間

conf.int = TRUE) # 增加置信區間

ggsurvplot(fit, data = lung,

conf.int = TRUE, # 增加置信區間

fun = "cumhaz") # 繪製累計風險曲線

ggsurvplot(fit, data = lung,

conf.int = TRUE, # 增加置信區間

risk.table = TRUE) #不同時間點風險人數表

ggsurvplot(fit, # 創建的擬合對象

data = lung, # 指定變量數據來源

conf.int = TRUE, # 顯示置信區間

pval = TRUE, # 添加P值

surv.median.line = "hv", # 添加中位生存時間線

add.all = TRUE) # 添加總患者生存曲線

六、採用log-rank 檢驗生存率差異

survdiff(Surv(time,status) ~ sex)

結果顯示：P=0.001,提示不同性別組生存率差異有統計學意義。

七、技能拓展

比較三組及三組以上生存時間和生存率的差異

7.1首先，查看分組變量ph.ecog,共四組，但有一組個數為1，故刪除該觀測，生成新的lung數據集。

table(ph.ecog)

lung<-lung[ph.ecog<3,]

7.2繪製曲線

fit <- survfit(Surv(time,status) ~ ph.ecog, # 創建生存對象

data = lung) # 數據集來源

ggsurvplot(fit, # 創建的擬合對象

data = lung, # 指定變量數據來源

conf.int = TRUE, # 顯示置信區間

pval = TRUE, # 添加P值

surv.median.line = "hv") # 添加中位生存時間線

7.3整體比較

survdiff(Surv(time,status) ~ ph.ecog,lung)

7.4兩兩比較

pairwise_survdiff(Surv(time,status) ~ ph.ecog,lung,p.adjust.method = "BH")

參考資料（可作為拓展閱讀材料）：

1. http://blog.sciencenet.cn/home.php?mod=space&uid=2609994&do=blog&id=976898

2. http://www.bio-info-trainee.com/1313.html

3. https://www.jianshu.com/p/3527d7bd6a49

4. https://www.cnblogs.com/wwxbi/p/6136348.html

5. https://zhuanlan.zhihu.com/p/49482538

6. https://blog.csdn.net/zfcjhdq/article/details/83502854

7. https://www.sohu.com/a/230205991_313188

8. https://www.jianshu.com/p/1e8b5286f007

9. R語言統計與繪圖：ggsurvplot()函數繪製Kaplan-Meier生存曲線

10. 王媛媛,畢玉,王在翔,宋棠,郭曉雷,付振濤,吳炳義. 山東省肺癌患者生存分析[J]. 中國衛生統計,2018,01:111-113+116.

統計諮詢

相關焦點

  • 生存曲線的 log-rank 檢驗
    生存曲線的比較：兩條或多條生存曲線的差別比較是生存分析的主要內容之一。
  • R數據分析:生存分析的做法和結果解釋
    也可以看到兩組隨時間變化的死亡比例是有顯著差異的，接下來寫寫不同生存曲線比較的檢驗：生存曲線的比較上面的例子中，我們分男女做了兩個生存曲線，這兩個生存曲線有沒有統計學差異呢？這時候就要用到log-rank test了：surv_diff <- survdiff(Surv(time, status) ~ sex, data = lung)surv_diff通過比較，我們發現兩個生存曲線確實存在顯著差異，此時我們就可以說性別為2的病人確實比性別為1的病人活得久點。
  • r 秩和檢驗 - CSDN
    所述配對雙樣品的Wilcoxon檢驗一種的非參數檢驗，其可以被用於比較樣品的兩個獨立數據。 本文介紹如何在ř中計算兩個樣本的秩檢驗。
  • 何使用Survminer包優雅的繪製生存曲線?
    引言：   生存分析是臨床上較為常用的統計學方法，用於比較不同組別的患者在接受幹預之後，生存時間的變化情況。生存分析是醫學領域中一個重要的內容，在各個疾病領域的研究中都運用十分廣泛。在R中進行生存分析常用的包主要有survival包以及survminer包。?
  • Stata教程第21期:K-M曲線繪製和Log-Rank檢驗
    上節課程我們一起學習了包括線性回歸、Logistic回歸、泊松回歸、Log-Binomial回歸在內的「廣義線性模型」，不知道小夥伴們是不是都完全消化了~ 本節課程，我們將通過4步，來詳解很多人都關注的「K-M曲線和
  • r語言 做wald檢驗_r語言wald檢驗怎麼做 - CSDN
    我們在做回歸擬合數據時，經常對因變量和自變量的假定是：自變量和因變量呈線性關係，logistic回歸是自變量和logitP，cox是指自變量和一個h(t)變換的一個東東（記不清了，抱歉），這些假定就叫線性假定。
  • r語言的p值檢驗 - CSDN
    )醫學統計R語言：分面畫boxplot醫學統計與R語言：調節效應分析(Moderation Analysis)醫學統計與R語言：結構方程模型(structural equation model)醫學統計與R語言：中介效應分析(mediation effect analysis)醫學統計與R語言：生存曲線(survival curves
  • 基於R的生存資料預測模型構建與Nomogram繪製
    根據上述條件，通過一定的數學模型判斷每個變量的得分，比如年齡40歲，得分是10分，男性得分是4分……依次累計各個變量的得分，計算出總分，不同的總分值就對應3月、6月和1年的生存概率。複雜的Cox回歸公式變成了直觀的圖片，醫生可以很方便的計算每個患者的生存概率，從而告知患者一個相對準確的「算命」結果。模型預測是否準確，一致性好不好等概念我們將在後續文章中陸續介紹。2.
  • R語言中繪製常用函數曲線的方法詳解
    在本號前面的文章中介紹了使用plot函數繪製圖形的方法。本文將使用另外一種方法來繪製常見函數的曲線的方法，如冪函數曲線、指數函數曲線、三角函數曲線等。所用函數簡介本文要使用到的R函數是curve函數。，分別是：y = x2和y = x2 + 3x + 1，其圖像見上圖藍色曲線。
  • R語言 小wald檢驗_lm檢驗 wald檢驗 - CSDN
    我們在做回歸擬合數據時，經常對因變量和自變量的假定是：自變量和因變量呈線性關係，logistic回歸是自變量和logitP，cox是指自變量和一個h(t)變換的一個東東（記不清了，抱歉），這些假定就叫線性假定。
  • R語言arma模型診斷_arma模型實現模型r語言 - CSDN
    library(nlme)          #調用其中的gls函數library(fArma)        #進行擬合和檢驗【基本函數】數學函數abs，sqrt：絕對值，平方根 log, log10, log2 , exp：
  • GraphPad Prism7繪製生存曲線
    GraphPad Prism7繪製生存曲線可以算得上是最簡單的一種方式。但是對於初學者來說，還是需要了解最基本的操作方法。今天，小編和大家一起來看看如何操作。首先，打開GraphPad Prism7，選擇「Survival」,如下圖所示在上圖中，小編標註了「1」和「2」，其中「1」是指輸入自己準備的數據。「2」是指使用軟體的示例數據。這裡小編就偷一下懶，使用示例數據。
  • r語言卡方檢驗和似然比檢驗_r語言似然比檢驗代碼 - CSDN
    本文對應《R語言實戰》第9章：方差分析；第10章：功效分析 ====================================================================方差分析：回歸分析是通過量化的預測變量來預測量化的響應變量，而解釋變量裡含有名義型或有序型因子變量時
  • 聊個天就把生存分析給做了?
    前些天我和我的學徒們在生信技能樹分享了教程：人人都可以學會生存分析（學徒數據挖掘），提到根據公共資料庫（比如TCGA）的信息來檢查具體的某個或者某幾個基因的生存意義已經是非常簡單了，只需要很簡單的資料庫網頁工具認知，或者一點點代碼能力。
  • r語言 似然比檢驗_對數似然比檢驗的r語言實現 - CSDN
    沒有必要做對比，因為我們不是在做兩兩比較。 結果表輸出類似於Wald檢驗結果，列與我們前面觀察到的相同。為什麼在LRT檢驗中會出現倍數變化?對於使用似然比檢驗的分析，p值僅由完整模型公式和簡化模型公式之間的偏差之差決定。
  • r語言tseries - CSDN
    library(FinTS) #調用其中的自回歸檢驗函數library(fGarch) #GARCH模型library(nlme) #調用其中的gls函數library(fArma) #進行擬合和檢驗二、【基本函數】數學函數abs，sqrt：絕對值，平方根 log, log10, log2 , exp：對數與指數函數
  • 醫學統計與R語言:GiViTI Calibration Belt
    >醫學統計與R語言：生存曲線（survival curves）with  risk.table醫學統計與R語言：如何比較兩種診斷試驗的靈敏度和特異度？醫學統計與R語言：你知道nomogram的points和total points怎麼算嗎？醫學統計與R語言：Cleveland dot plot醫學統計與R語言：交互作用模型中分組效應及標準誤的計算醫學統計與R語言：多條ROC曲線的AUC多重比較醫學統計與R語言：來，今天學個散點圖！
  • r語言檢驗序列相關 - CSDN
    ，則進行模式識別畫自相關圖和非自相關圖，根據兩圖的結尾性和拖尾性進行AR、MA、ARMA的模式識別對識別後模式中的位置參數進行參數估計arima()模型檢驗分為：①殘差的白噪聲檢驗；②過度擬合檢驗pt()模型檢驗通過則進行模型優化，否則重新進行模式識別模型優化中得到AIC和BIC值，進行模型的優化然後進行預測與控制2.
  • R語言繪圖 | R繪製火山圖 EnhancedVolcano+ggplot
    一般默認取log2FC絕對值大於1為差異基因的篩選標準。Y軸：-Log(adjust P-value), 對矯正後的P值取負對數(-log)；矯正P值為多重假設檢驗矯正過的差異顯著性P值。由於轉錄組測序的差異表達分析是對大量的基因表達值進行獨立的統計假設檢驗會存在假陽性問題，因此在進行差異表達分析過程中，採用了公認的Benjamini-Hochberg校正方法對原有假設檢驗得到的顯著性p值（p-value）進行校正，剔除假陽性。