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

題記：本文成形於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.