基於R的生存資料預測模型構建與Nomogram繪製

2021-02-15 臨床研究與醫學統計

題記：本文成形於2015年並於2017年收錄於由筆者與胡志德醫生共同主編，中南大學出版社出版的《瘋狂統計學》，此次屬舊文重發，較原文有較大改動。

1. 背景介紹

人們對於「算命」這樣的事情總是樂此不彼，無論是中國文化中「算命佔卜」還是西方文化中的「佔星術」，無不說明人們對於此道的熱衷。本章討論的構建預測模型也是一種「算命」，只是這是一種更為科學的「算命」。作為一個腫瘤科醫生，臨床上可能會遇到這樣的情況，一個55歲男性中晚期食管癌患者問道：醫生，請問我還能活多久？大多數時候我們會根據患者的臨床分期告訴患者，這種情況的中位生存時間是多少……其實這裡的臨床分期就是我們預測患者期望壽命的依據，或者如本文標題所說的「預測模型」，我們用某一臨床分期的所有患者中位生存時間來回答單個患者的問題，可能會存在問題，這裡的中位生存時間代表的是患病人群預期壽命的平均值，對於具體的個體來說有時並不準確。

那我們是否有更科學、更準確的方法計算出單個患者的生存概率？答案是肯定的。我們可以通過Cox回歸等方法先構建一個數學模型，然後把這個模型可視化為Nomogram，通過圖形來計算每個具體患者的生存概率。Nomogram中文常稱為諾莫圖或者列線圖，其本質就是回歸方程的可視化。它根據所有自變量標準回歸係數的大小來制定評分標準，給每個自變量的每種取值水平一個評分，對每個患者，就可計算得到一個總分，再通過得分與結局發生概率之間的轉換函數來計算每個患者的結局時間發生的概率。舉例來說，我們治療了一個40歲的男性胰腺癌手術後的患者，術中進行了放療，腫瘤位於胰腺頭部，有腹膜轉移，臨床分期在IV期。根據上述條件，通過一定的數學模型判斷每個變量的得分，比如年齡40歲，得分是10分，男性得分是4分……依次累計各個變量的得分，計算出總分，不同的總分值就對應3月、6月和1年的生存概率。複雜的Cox回歸公式變成了直觀的圖片，醫生可以很方便的計算每個患者的生存概率，從而告知患者一個相對準確的「算命」結果。模型預測是否準確，一致性好不好等概念我們將在後續文章中陸續介紹。

2. 案例與操作

本章我們以[案例 1]的數據，介紹生存資料預後模型的建立與Nomogram圖形繪製。為了便於讀者閱讀與練習模仿，筆者對數據進行了簡化處理。


[案例 1] 筆者在The Cancer Genome Atlas (TCGA) 數據中下載1215例浸潤性乳腺癌患者的臨床資料。下載網址: https://genome-cancer.ucsc.edu/。數據經整理後如下 表1 所示，變量定義及賦值說明如表2所示。研究目的就是要根據1215例乳腺癌患者的生存資料構建預測模型並繪製 Nomogram。

表1. 1215例乳腺癌患者的生存資料

表2.變量定義、賦值及說明

整個預測模型構建的步驟大體如下：

步驟 1. 我們首先使用Cox回歸基於構建預測模型並篩選獨立預後因素（用於建模的數據集一般稱為訓練集或者內部數據集）。需要說明的是本例的數據錄入、單因素Cox回歸分析，多因素Cox回歸分析等操作可參考筆者主編的《聰明統計學》[1] 與《瘋狂統計學》[2]。最終我們可得到三個與預後相關的獨立因素：Age, PgR, Pathologic_stage。

步驟 2. 我們就以這三個獨立的預後因素繪製Nomogram，建模完成。

步驟 3. 對上述兩步所構建的預測模型的區分能力 (Discrimination) 進行評價，並計算C-index。

步驟 4. 對模型進行驗證，可通過外部數據集進行驗證，如果無法獲得外部數據集，筆者推薦採用Bootstrap衝抽樣法基於訓練集驗證模型並繪製校正曲線(Calibration plot)。

3. 計算過程

我們重點講解Cox回歸模型中C-Index計算，Nomogram繪製過程，Bootstrap法驗證模型及繪製標準曲線。以下所有操作均基於R軟體，軟體下載地址：https://www.r-project.org/. 我們把數據另存為BreastCancer.sav，並保存至R的當前工作目錄下。計算結果如圖1及圖2所示。

載入程序包

library(foreign)library(rms)

## Loading required package: Hmisc

## Loading required package: lattice

## Loading required package: survival

## Loading required package: Formula

## Loading required package: ggplot2

##
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:base':
##
##     format.pval, units

## Loading required package: SparseM

##
## Attaching package: 'SparseM'

## The following object is masked from 'package:base':
##
##     backsolve

數據準備，讀入當前工作路徑下的.sav格式的數據

breast<-read.spss("BreastCancer.sav")

把數據轉換為數據框格式

breast<-as.data.frame(breast)
breast<-na.omit(breast)

展示數據框的前6行數據

head(breast)

##    No    Months Status Age       ER      PgR Margin_status
## 9   9  8.633333 Censor  70 Positive Negative      Nagative
## 11 11 44.033333 Censor  56 Positive Positive      Nagative
## 12 12 48.766667 Censor  54 Positive Positive      Nagative
## 13 13 14.466667 Censor  61 Positive Positive      Nagative
## 14 14 47.900000 Censor  39 Negative Positive      Nagative
## 19 19 39.866667 Censor  50 Positive Positive      Nagative
##    Pathologic_stage HER2_Status Menopause_status
## 9           Stage I    Negative   Post menopause
## 11          Stage I    Negative    Pre menopause
## 12         Stage II    Negative    Pre menopause
## 13         Stage II    Negative   Post menopause
## 14         Stage II    Negative    Pre menopause
## 19         Stage II    Positive   Post menopause
##                 Surgery_method             Histological_type
## 9                   Lumpectomy                         Other
## 11 Modified Radical Mastectomy                         Other
## 12 Modified Radical Mastectomy Infiltrating Ductal Carcinoma
## 13                  Lumpectomy                         Other
## 14                  Lumpectomy Infiltrating Ductal Carcinoma
## 19                  Lumpectomy Infiltrating Ductal Carcinoma

定義終點事件

breast$Status<-ifelse(breast$Status=="Dead",1,0)

設置參考組

breast$Pathologic_stage<- relevel(breast$Pathologic_stage,ref='Stage I')

構建Cox回歸方程

coxm <-cph(Surv(Months,Status==1)~Age+Pathologic_stage+PgR,x=T,y=T,data=breast,surv=T)
surv<- Survival(coxm) surv1<- function(x)surv(1*12,lp=x) surv2<- function(x)surv(1*36,lp=x) surv3<- function(x)surv(1*60,lp=x) dd<-datadist(breast) options(datadist='dd')

繪製列線圖。maxscale參數指定最高分數，一般設置為100或者10分；un.at 設置生存率的刻度；xfrac 設置數值軸與最左邊標籤的距離，可以調節下數值觀察下圖片變化情況

plot(nomogram(coxm,fun=list(surv1,surv2,surv3),lp=F,funlabel=c("1-Yeas OS", '3-Year OS','5-YearOS'),maxscale=100,fun.at=c('0.95','0.85','0.80','0.70','0.6','0.5','0.4','0.3','0.2','0.1')),xfrac=.3)

圖1. 列線圖。

結果解讀：圖中points就是一個選定的評分標準或者尺度，對於每個自變量取值，在該點做一條垂直於Points軸的直線（可通過直尺），交點即代表該自變量取值下的評分，如Age，取25時評分為0分，CEA為90時評分約為100分。以此類推，計算每個患者各個自變量對應的分值points，加起來就是總分totalpoints。同樣的道理，在TotalPoints軸上找到該患者總分對應的點，畫一垂直直接到生存概率軸上（如3年生存概率3-yearOS或者5年生存概率5-year OS），交點即為該患者的3年或5年生存概率。

計算C-index

f<-coxph(Surv(Months,Status==1)~Age+Pathologic_stage+PgR,data=breast)
sum.surv<-summary(f)
c_index<-sum.surv$concordance
c_index

##          C      se(C)
## 0.77878820 0.05734042

此處再次解釋一下上述R語言代碼中計算的C-Index統計量的含義，此值可以參考ROC曲線下面積去理解，取值範圍0~1，越接近1說明我們構建的Cox回歸模型預測價值越大，一般來說C-Index=0.7即表明模型具有良好的預測價值。本例中計算的C-Index = 0.7503, se(C-Index) = 0.02992，以上為軟體直接輸出結果。

計算C-index的補數（可選）

library(Hmisc)
S<-Surv(breast$Months,breast$Status==1)
rcorrcens(S~predict(coxm),outx=TRUE)

##
## Somers' Rank Correlation for Censored Data    Response variable:S
##
##                  C    Dxy  aDxy   SD    Z P   n
## predict(coxm) 0.22 -0.559 0.559 0.09 6.21 0 549

構建標準曲線。u需要與之前f中定義好的time.inc一致，即f中time.inc為60時，此處u即為60；m要根據樣本量來確定，由於標準曲線一般將所有樣本分為3或4組（在圖中顯示3或4個點），m代表每組的樣本量數，因此m*3應該等於或近似等於樣本量。

cal<- calibrate(coxm, cmethod='KM', method='boot', u=60, m=100, B=100)

## Using Cox survival estimates at  30 Days

## Warning in groupkm(cox, Surv(y[, 1], y[, 2]), u = u, cuts = orig.cuts): one
## interval had < 2 observations

列印並修改標準曲線的圖形參數

plot(cal,lwd=2,lty=1,errbar.col=c(rgb(0,118,192,maxColorValue=255)),xlim=c(0.7,1),ylim=c(0.6,1),xlab="Nomogram-Predicted Probability of 5-Year OS",ylab="Actual 5-Year OS(proportion)", col=c(rgb(192,98,83,maxColorValue=255)))
lines(cal[,c("mean.predicted","KM")],type="b",lwd=2,col=c(rgb(192,98,83,maxColorValue=255)), pch=16)
abline(0,1,lty=3,lwd=2,col=c(rgb(0,118,192,maxColorValue=255)))

圖2. 標準曲線。

基於訓練集採用bootstrap重抽樣的方法驗證列線圖的預測準確性。結果解讀：橫坐標為根據列線圖預測的每個患者的5年生存概率，縱坐標為每個患者實際的5年生存概率，如果圖中的紅色線條恰好與藍色虛線完全重合時最為理想。

4. 總結與討論

綜上，本文介紹了生存資料預測模型的構建以及Nomogram的繪製。需要說明的是一個預測模型是否具有較好的實際運用價值，除了考量其預測的準確性還要考量其可操作性。除了內部驗證檢驗其準確性，外部驗證有時也是必須的，本例中因未獲得較好的外部驗證數據及並未演示外部驗證的操作過程。過去十年有關臨床預測模型的Nomogram文章很多，很多Nomogram的預測效能遠高於TNM分期，可是為什麼醫生們仍然習慣運用TNM分期去為患者「算命」呢？我想主要原因還是TNM分期使用更為方便。從這個角度講，我們構建Nomogram應該納入儘可能少的變量以方便臨床使用。到底應該優先考慮其準確性還是應該考慮其可操作性？讀者應該根據自己的研究目的自行考量。

5. 參考文獻

[1]. 周支瑞, 胡志德. 聰明統計學. 長沙:中南大學出版社, 2016.

[2]. 周支瑞, 胡志德. 瘋狂統計學. 長沙:中南大學出版社, 2018.

相關焦點

  • 天津醫科大學總醫院王邦茂、孫超團隊:基於腰大肌指數構建肝硬化患者長期預後的nomogram預測模型 | AME作者面對面
    所刊雜誌：Annals of Translational Medicine（點擊查看雜誌詳情與影響因子）文章標題：基於腰大肌指數構建肝硬化患者長期預後的nomogram預測模型（Low psoas muscle index associates with
  • 教你用R畫列線圖(Nomogram),讓預測模型結果可視化!
    還沒來得及閱讀的小夥伴請點擊查看：　　　　同樣是構建多因素回歸模型，往往我們另一個主要目的是為了對結局事件的發生風險進行預測，那麼是否也可以將預測模型的結果，像森林圖那樣可視化地展示出來呢？今天小咖就來帶大家認識一下神奇的列線圖。
  • 瘋狂統計學 | 生存資料預後模型的建立及Nomogram的繪製方法
    事實上本章討論的構建預後模型也是一種「算命」，只是這是一種更為科學的 「算命」。筆者作為一個腫瘤專科醫生，在日常工作中也經常會遇到如下的情況：一位55歲的男性中晚期食管癌患者絕望地詢問我：「周醫生，請問我還能活多久？」
  • r語言 做wald檢驗_r語言wald檢驗怎麼做 - CSDN
    instead of x1:x2 because x1 %ia% x2 triggers# anova to pool x1*x2 term into x1 terms to test total effect# of x1## End(Not run) SVM的調參，關於e1071包，好像如果把數據「尺度化」（scaling）後，使用默認的參數就能訓練出較好的模型
  • 耿智敏|膽囊癌相關生存預測模型的研究現狀
    [7]羅勝蘭.生存分析的方法及應用[J].浙江預防醫學, 2013, 25(5):29-31,34.DOI:10.3969/j.issn.1007-0931.2013.05.008.[8]孫振球，徐勇勇.醫學統計學[M].4版.北京:人民衛生出版社,2014:291-294.
  • R語言 小wald檢驗_lm檢驗 wald檢驗 - CSDN
    instead of x1:x2 because x1 %ia% x2 triggers# anova to pool x1*x2 term into x1 terms to test total effect# of x1## End(Not run) SVM的調參，關於e1071包，好像如果把數據「尺度化」（scaling）後，使用默認的參數就能訓練出較好的模型
  • 上海交大團隊提出基於信號通路的肝癌風險預測模型
    近日，上海交通大學生命科學技術學院、生物信息學與生物統計學系俞章盛教授團隊在肝癌風險預測研究方面取得新進展，在柳葉刀雜誌子刊《EBioMedicine》發表了題為「Pathway-based Biomarker Identification with Crosstalk Analysis for Robust Prognosis Prediction
  • R語言arma模型診斷_arma模型實現模型r語言 - CSDN
    acf(sha,22)   #繪製自相關圖，滯後期數22pacf(sha,22)  #繪製偏自相關圖，滯後期數22corr=acf(sha,22)   #保存相關係數cov=acf(sha,22,type = "covariance")   #保存協方差
  • 你為何而轉:微博用戶轉發行為預測模型的構建與影響因素探究
    本研究旨在構建模型對轉發行為進行預測，並分析其影響因素。首先根據「5W模型」，提取微博作者、微博文本、微博接受者和相互關係四個一級維度特徵，並細分為39個二級維度特徵，利用支持向量機構建預測模型，再通過新浪微博數據，對模型進行評估。預測模型的查全率為58.67%，精確率為82.19%，F1值為68.46%，這表明預測模型的表現令人滿意。
  • 用ROC曲線解析預測模型
    最具吸引力的目標之一是利用數據資產的力量來創建能夠預測各種結果的機器學習模型。通過定義明確的模型，可以確定能夠預測結果的最重要影響因素，為戰略假設開發有價值的洞察力，甚至可以通過友好的用戶界面將模型的邏輯實現到軟體應用程式中。然而，在任何這種魔法發生之前，我們需要知道模型創建的預測是否有益！
  • 通過Keras 構建基於 LSTM 模型的故事生成器
    所以神經網絡要準確進行預測，就必須記憶單詞的所以序列。 而這正是 LSTM 可以做到的。編程實現 LSTM本文將通過 LSTM 網絡開發一個故事生成器模型。主要使用自然語言處理（NLP）進行數據預處理，使用雙向LSTM進行模型構建。
  • 何使用Survminer包優雅的繪製生存曲線?
    Survival 包提供了生存函數的建立，Cox模型的建立，以及比較分析。這個包也提供了基於基礎繪圖系統的生存曲線繪製。?Survminer包提供了基於ggplot2系統對於生存分析的可視化，使得生存分析具有更加美觀的圖形，以及自我定製方式。而今天我們的主角就是Survminer包，讓我們鼓足精神一起來學習「如何使用Survminer包優雅的繪製生存曲線」吧。
  • 第四十講 R-線性回歸:預測模型及可信區間
    今天的課程將繼續帶大家學習R-線性回歸：預測模型及可信區間。線性回歸的一個主要目標是基於一個或多個預測變量來預測結果值。（我們也用它來研究兩個變量的相關性，同時校正其他混雜因素）。那麼，當我們取得了預測模型後，根據該預測模型對新數據進行預測得出的預測值是什麼？這個預測值的可信度如何呢？今天的講解中，我們會給出答案。我們首先建立一個簡單的線性回歸模型，該模型根據速度預測汽車的停車距離。
  • 基於DNA甲基化的分子亞型構建發5+分
    ，並根據亞型特異甲基化位點構建COAD的COX比例風險模型預測患者預後，快來學習一下吧!作者想要通過對COAD組織甲基化晶片的分析，根據樣本的甲基化水平將其區分為不同的分子亞型；並且基於COAD患者癌組織甲基化水平建立預後預測模型，以改善對COAD預後的評估。二. 文章思路三.
  • 財務風險預警模型構建實證分析
    【摘要】財務風險預警模型可以從定量角度客觀準確判斷企業的財務危機程度，網絡的普及對此頗具影響。本文採用實證方法，選取了5個財務指標作為模型變量，構建了一個基於極值原理的Fisher線性判別模型，並對該模型進行了實證檢驗。結果表明，該模型具有較好的對企業財務狀況和風險狀況進行評價預警的能力。
  • 時間序列的LSTM模型預測——基於Keras
    時間序列的預測指的是基於序列的歷史數據，以及可能對結果產生影響的其他相關序列，對序列未來的可能取值做出預測。現實生活中的時間序列數據預測問題有很多，包括語音分析、噪聲消除以及股票期貨市場的分析等，其本質主要是根據前Ｔ個時刻的觀測數據推算出Ｔ＋１時刻的時間序列的值。
  • 不用代碼,教你Excel構建數據分析預測模型!
    下面是一個教程，介紹如何在Excel中構建線性回歸模型以及如何解釋結果。Excel真的能構建預測模型？這通常是我提起這個話題時的第一反應。當我演示如何利用Excel的靈活性為我們的數據科學和分析項目構建預測模型時，接下來是一個令人懷疑的眼神。讓我問你一個問題：如果你周圍的商店開始收集客戶數據，他們是否可以採用基於數據的策略來銷售他們的商品？
  • 從這篇22分+文章入手,帶你深度探討臨床預測模型研究思路
    臨床預測模型是基於疾病的特徵表型建立統計模型，用來預測具有某些特徵的人群未來某種結局事件發生的概率。主要包括疾病診斷預測模型與疾病預後預測模型兩種形式。預測模型構建思路選擇建模方法預測模型構建的方法主要包括參數化模型，半參數化模型和非參數化模型。
  • 中國大陸基于震源物理模型的長期地震預測及其展望
    報告題目：中國大陸基于震源物理模型的長期地震預測及其展望報 告 人：邵志剛 研究員報告時間：2020年12月4日（周五）上午9：00-11：30報告地點：地質宮519報告摘要：以中國大陸活動地塊理論為指導，長期地點預測重點關注活動地塊邊界帶，綜合強震破裂空段、斷層運動閉鎖段
  • 基於GWO-BP-CNN-ec的風電功率短期預測模型*
    通過仿真分析得出，基於GWO-BP-CNN-ec的組合模型預測效果優於單一CNN，GWO-CNN及GWO-BP-CNN，然後對比了基於BP、Elman、小波神經網絡的淺層網絡預測模型與基於CNN的深層網絡組合預測模型的預測效果，驗證出了深層CNN網絡在預測方面的精度優勢，在高風速段時尤為明顯。證明深層卷積神經網絡組合模型對風電短期功率預測精度的提升有顯著作用。