【原】生存資料（參數(shù)、半?yún)?shù)）預(yù)測模型的列線圖繪制

Memo_Cleon 2023-03-06 發(fā)布于上海

展開全文

列線圖（Alignment Diagram），也稱諾莫圖(Nomogram），絕是臨床預(yù)測模型最為常見的一種呈現(xiàn)形式。列線圖可以根據(jù)受試者的一些特征（解釋變量）來預(yù)測結(jié)局事件的發(fā)生概率或者患者的生存率等。生存資料預(yù)測模型的列線圖除了可以預(yù)測生存率，還可以預(yù)測分位數(shù)生存時(shí)間、生存時(shí)間期望均值、風(fēng)險(xiǎn)比。

生存資料預(yù)測模型的構(gòu)建絕大多數(shù)是基于COX回歸法。COX回歸屬于半?yún)?shù)法，在多因素的生存分析中幾乎是霸主般的存在，主要是因?yàn)?/span>COX回歸幾乎不需要考慮前提條件，僅要求滿足比例風(fēng)險(xiǎn)假定就可以了。除了COX回歸，參數(shù)生存模型也是生存資料非常重要的分析方法。但參數(shù)法需要數(shù)據(jù)滿足一定分布，一說到分布啥的一些理論，基本上就阻斷了大部分人繼續(xù)下去的欲望，甚至頭皮發(fā)麻惡心想吐，哪還有心思去搞清楚怎么用。但其實(shí)如果滿足適用條件，參數(shù)法可能會(huì)比非參數(shù)法更好地解釋數(shù)據(jù)。正如在比較組間效應(yīng)時(shí)，如果數(shù)據(jù)滿足獨(dú)立正態(tài)同方差，我們更愿意適用t檢驗(yàn)/方差分析，而不會(huì)首選非參數(shù)檢驗(yàn)一樣，雖然使用非參數(shù)檢驗(yàn)也沒有錯(cuò)。與生存時(shí)間相關(guān)的一些重要分布有"weibull", "exponential", "gaussian", "logistic","lognormal" and "loglogistic"。

此次筆記的任務(wù)就是繪制生存資料預(yù)測模型的列線圖，包括COX回歸模型和參數(shù)生存模型。

數(shù)據(jù)來源：UCI Machine Learning Repository

[https://archive.ics./ml/datasets/Heart+failure+clinical+records#]

這是一個(gè)心力衰竭臨床記錄數(shù)據(jù)集，有13個(gè)臨床特征。筆者對(duì)個(gè)別字段的數(shù)據(jù)類型進(jìn)行了轉(zhuǎn)換。

下面直接上命令清單：

導(dǎo)入數(shù)據(jù)與因子設(shè)置

##導(dǎo)入數(shù)據(jù)

library(haven)

HF<- read_sav("D:/Temp/HF.sav")

HF<-as.data.frame(HF)

##設(shè)置因子及因子各水平的標(biāo)簽

HF$Age.cat<-factor(HF$Age.cat,levels=c(0,1,2),labels=c("<60","60-80",">=80"))

HF$SCr.cat<-factor(HF$SCr.cat,levels=c(0,1),labels=c("Normal","Higher"))

……

【1】COX回歸模型列線圖的繪制

1.1 建立COX回歸模型

在經(jīng)過單變量分析和全子集回歸的分析后（此處略），最終選擇了三個(gè)變量（Age.cat、EF、SCr.cat）來建立COX生存模型，模型正態(tài)基本滿足PH假定。

1.2.1 COX回歸模型列線圖繪制

##載入列線圖繪制程序包

library(rms)

library(Hmisc);library(lattice);library(survival);library(Formula);library(ggplot2);library(SparseM)

##設(shè)置工作環(huán)境，匯總預(yù)測變量分布信息。該設(shè)置為必需項(xiàng)。

CoxNM<-datadist(HF)

options(datadist='CoxNM')

##設(shè)置預(yù)測變量在列線圖上顯示的名稱。不設(shè)置則默認(rèn)為變量的名稱

label(HF$Age.cat)<-"Age"

label(HF$EF)<-"Ejection Fraction"

label(HF$SCr.cat)<-"Serum Creatinine"

##建立COX回歸模型

HFDeath.cox<-cph(Surv(Time,Death==1)~Age.cat+EF+SCr.cat,data=HF,x=T,y=T,surv=T,method="breslow")

#method：結(jié)（ties）的處理方式有efron、breslow和exact，默認(rèn)為efron。這里改為method="breslow"是為了得到跟SPSS一致的結(jié)果，STATA默認(rèn)的也是breslow。同樣的數(shù)據(jù)在采用R的coxph{survival}和cph{rms}的進(jìn)行Cox回歸，其結(jié)果與SPSS有略微差異，如果你像我一樣糾結(jié)，可以在R中把結(jié)的處理方式改一下就可以了；

#x、y、surv：這幾個(gè)參數(shù)是為了繪制后面的列線圖方便，都設(shè)置為TRUE就可以了；

#time.inc：這在后期進(jìn)行模型評(píng)估中可能會(huì)用到，比如評(píng)估模型在不同時(shí)間點(diǎn)上的校準(zhǔn)度，常常用函數(shù)calibrate繪制1年、3、5年的校準(zhǔn)曲線，這時(shí)候就需要用time.inc指定時(shí)間點(diǎn)，然后再calibrate函數(shù)的u參數(shù)進(jìn)行呼應(yīng)。

##建立生存率的估計(jì)函數(shù)

surv<-Survival(HFDeath.cox)

surv3M<-function(x)surv(91.31,lp=x)

surv0.5<-function(x)surv(182.625,lp=x)

surv1<-function(x)surv(285,lp=x)

#Survival函數(shù)可以對(duì)cph創(chuàng)建的對(duì)象返回一個(gè)S函數(shù)，用于計(jì)算生存率的估計(jì)值。除了Survival函數(shù)，Quantile()返回的S函數(shù)可用于計(jì)算生存時(shí)間的分位數(shù)（默認(rèn)為中位數(shù)，可通過q來指定具體的分位數(shù)）；Mean()返回返回一個(gè)計(jì)算平均生存時(shí)間的函數(shù)。參數(shù)估計(jì)模型還有Hazard()用于估計(jì)風(fēng)險(xiǎn)函數(shù)；

#通過function函數(shù)定義了三個(gè)生存率的計(jì)算函數(shù)，這些函數(shù)以一個(gè)時(shí)間點(diǎn)和線性預(yù)測器作為參數(shù)，并返回該時(shí)間點(diǎn)的生存率。本例時(shí)間單位是天，用91.31天和182.62天來表示1個(gè)月和3個(gè)月，1年應(yīng)該用365.25天，但本例最長的隨訪時(shí)間是285天，如果使用大于285的數(shù)字，結(jié)果會(huì)報(bào)錯(cuò)“Error in approx(fu[s], xseq[s], fat, ties = mean) :調(diào)用內(nèi)插至少需要兩個(gè)非NA值的數(shù)”。此處僅示例。

##列線圖的參數(shù)設(shè)置

Nom.Surv<-nomogram(HFDeath.cox,lp=T,fun=list(surv3M,surv0.5,surv1), fun.at=c('0.99','0.95','0.90','0.85',seq(0.80,0.10,by=-0.1),'0.05','0.01'),funlabel=c("3-Month Survival Prob.","Half-Year Survival Prob.","1-Year Survival Prob."))

##列線圖繪制的主要函數(shù)，參數(shù)又多又復(fù)雜又難理解，僅介紹幾個(gè)常用的、重要的

#lp：默認(rèn)設(shè)置為FALSE，即不顯示線性預(yù)測器（linear predictor）的那條坐標(biāo)軸。我們最終想預(yù)測的概率就是這個(gè)線性預(yù)測器的值通過一定的函數(shù)變換得來的。列線圖的基本原理就是根據(jù)預(yù)測變量對(duì)結(jié)局的影響大?。ɑ貧w系數(shù)），給各個(gè)預(yù)測變量的水平進(jìn)行賦分，評(píng)分之和通過函數(shù)轉(zhuǎn)化成結(jié)局事件發(fā)生概率。這個(gè)線性預(yù)測值可以理解為承接變量總評(píng)分與預(yù)測概率之間的中間變量，對(duì)列線圖的使用而言其實(shí)沒啥毛用；

#fun：對(duì)線性預(yù)測值進(jìn)行轉(zhuǎn)換的函數(shù)；

#fun.at：預(yù)測值橫軸的刻度，本例是在坐標(biāo)軸上顯示的概率值?？墒褂迷搮?shù)的默認(rèn)值，也可通過該參數(shù)改成需要顯示的值。同一個(gè)fun函數(shù)可以是使用同一個(gè)fun.at向量值，不同fun函數(shù)需要用fun.at值的向量列表標(biāo)識(shí)。本例有三個(gè)自建函數(shù)，不過都是對(duì)應(yīng)生成生存率的Survival這一個(gè)函數(shù)，因此三條預(yù)測概率坐標(biāo)軸刻度使用了相同的fun.at值;

#funlabel：預(yù)測橫軸名稱，默認(rèn)"Predicted Value"；

#conf.int: 是否顯示得分的置信區(qū)間，默認(rèn)否；

#maxscale：默認(rèn)單變量最大分?jǐn)?shù)值是100分。

##繪制、美化列線圖

plot(Nom.Surv,xfrac=.25,col.grid=gray(c(0.8,0.95)))

##主要是將前面Nomogram函數(shù)的結(jié)果畫出來，一些參數(shù)可以起到美化作用

#lplabel：線性預(yù)測器橫軸的名稱，默認(rèn)"Linear Predictor"。不過這個(gè)變量我一般會(huì)在Nomogram中關(guān)掉不顯示；

#fun.side：預(yù)測軸上的刻度值位置，1用于定位在軸下方（默認(rèn)值），3用于軸上方。預(yù)測橫軸的刻度值有時(shí)候會(huì)非常擁擠甚至重疊在一起，如一個(gè)軸有5個(gè)標(biāo)記標(biāo)簽值，第二個(gè)和第三個(gè)相互重疊，可該參數(shù)進(jìn)行指定：fun.side=c(1,1,3,1,1)；

#col.conf：置信區(qū)間的顏色設(shè)置；

#xfrac：各刻度線橫軸與軸名稱之間的空白間隔所占的比例，默認(rèn)0.35。如果你覺得坐標(biāo)軸名稱與坐標(biāo)軸隔得距離太遠(yuǎn)，可以把這個(gè)值設(shè)置的小一些；

#cex.axis：刻度軸刻度標(biāo)簽的字符大小，默認(rèn)0.85；

#cex.var：刻度軸名稱的字符大小，默認(rèn)1；刻度軸標(biāo)題默認(rèn)是變量名稱，可通過label()進(jìn)行修改；

#col.grid ：默認(rèn)不繪制垂直參考網(wǎng)格線。該參數(shù)若長度為1，則對(duì)主要和次要參考線使用相同的顏色，長度為2則分別對(duì)應(yīng)主要和次要參考線的顏色，常用灰色col.grid = gray(c(0.8, 0.95))；

#points.label：單變量評(píng)分軸的名稱，默認(rèn)為'Points'；

#total.points.label：總評(píng)分軸的名稱，默認(rèn)'Total Points'。

列線圖的主要控制參數(shù)及繪制結(jié)果如下，列線圖的解讀方法可參見第2部分參數(shù)生存模型的列線圖部分。

1.2.2帶有交互作用的列線圖繪制

模型是否存在交互作用應(yīng)在1.2步進(jìn)行考察，此處僅演示其列線圖的繪制操作。

##上接設(shè)置工作環(huán)境變量顯示名稱

HFDeath.cox2<-cph(Surv(Time,Death==1)~Age.cat+EF+SCr.cat+Age.cat*EF,data=HF,x=T,y=T,surv=T,method="breslow")

surv<-Survival(HFDeath.cox2)

mea<-Mean(HFDeath.cox2)

surv3M<-function(x)surv(91.31,lp=x)

surv0.5<-function(x)surv(182.625,lp=x)

surv1<-function(x)surv(285,lp=x)

mean<-function(x)mea(lp=x)

#除了建立生存率預(yù)測函數(shù)，還建立一個(gè)生存時(shí)間期望均值的預(yù)測函數(shù)。最后一次隨訪時(shí)結(jié)局為非刪失事件，Mean函數(shù)會(huì)有報(bào)警。由于刪失的存在，平均生存時(shí)間是一個(gè)有偏估計(jì)量，一般用生存曲線下面積或累計(jì)發(fā)病率曲線上面積表示。

Nom.Surv2<-nomogram(HFDeath.cox2,lp=F,fun=list(surv3M,surv0.5,surv1,mean),fun.at=c(c(seq(0.90,0.10,by=-0.1),'0.05'),seq(280,20,by=-20)),funlabel=c("3-Month Survival Prob.","Half-Year Survival Prob.","1-Year Survival Prob.","Mean Survival Time"))

plot(Nom.Surv2,xfrac=.25,col.grid=c(7,3))

#此處本來還想演示一下fun.side參數(shù)的：fun.side=list(rep(1,10),rep(3,10),rep(c(1,3),5),c(rep(1,10),3,1,3,1))，結(jié)果報(bào)錯(cuò)。單獨(dú)預(yù)測生存率或者單獨(dú)預(yù)測生存時(shí)間均值時(shí)都沒有問題，但合在一起就出現(xiàn)了錯(cuò)誤，尋因良久，未果，只得作罷。不過可以適當(dāng)對(duì)fun.at參數(shù)進(jìn)行調(diào)整來解決刻度值擁擠的問題。

對(duì)某一個(gè)患者，在預(yù)測時(shí)交互作用只選擇其中符合該患者特征那一行進(jìn)行評(píng)分即可，比如要預(yù)測的患者70歲，EF和Age.cat的交互作用我們只用列線圖EF(Age.cat=60-80)這一個(gè)橫軸，結(jié)合EF值來打分就可以了。列線圖的使用可參見第2部分參數(shù)生存模型的列線圖部分。

【2】參數(shù)生存模型

rms包中的psm函數(shù)可以提供6種參數(shù)生存模型（Parametric Survival Model）的擬合，包括威布爾分布（weibull），指數(shù)分布（exponential）、正態(tài)分布（gaussian）、logistic、對(duì)數(shù)正態(tài)（lognormal）和對(duì)數(shù)logistic。

經(jīng)初步考察，在psm函數(shù)中提供的6個(gè)參數(shù)分布中，數(shù)據(jù)擬合對(duì)數(shù)正態(tài)分布似乎更好一些。

##上接設(shè)置工作環(huán)境變量顯示名稱

HFDeath.psm<-psm(Surv(Time,Death==1)~Age.cat+EF+SCr.cat,data=HF, dist="lognormal")

#dist：可選擇的分布有weibull, exponential、gaussian、logistic、lognormal和loglogistic

surv<-Survival(HFDeath.psm)

quan<-Quantile(HFDeath.psm)

surv100D<-function(x)surv(100,lp=x)

surv0.5<-function(x)surv(182.625,lp=x)

surv1<-function(x)surv(285,lp=x)

Q25<-function(x)quan(lp=x,q=0.25)

median<-function(x)quan(lp=x)

#除了建立生存率預(yù)測函數(shù)，同時(shí)建立生存時(shí)間分位數(shù)的預(yù)測函數(shù)。默認(rèn)是中位值，通過q參數(shù)可指定任意分位數(shù)

#本例25分位數(shù)生存時(shí)間為100天，隨訪結(jié)束時(shí)的累積生存率為57.6%，尚未到達(dá)中位生存時(shí)間

Nom.Surv.psm<-nomogram(HFDeath.psm,lp=F,fun=list(surv100D,surv0.5,surv1,Q25,median), fun.at=c(c('0.05',seq(0.10,0.90,by=0.1)),c(seq(20,100,by=20),seq(150,300,by=50),seq(400,1200,by=200),2000,2500)),funlabel=c("100-Days Survival Prob.","Half-Year Survival Prob.","1-Year Survival Prob.","Q25 Survival Time","Q50 Survival Time"))

plot(Nom.Surv.psm,xfrac=.25,col.grid=c("red","green"))

#采用默認(rèn)的fun.at，分位數(shù)預(yù)測刻度軸的標(biāo)簽會(huì)非常得擁擠，但在調(diào)整展示兩個(gè)函數(shù)（此處是生存函數(shù)和分為數(shù)函數(shù)）刻度值的位置時(shí)，fun.side參數(shù)再次報(bào)錯(cuò)，只能使用fun.at可以進(jìn)行適當(dāng)?shù)恼{(diào)整。

一個(gè)年齡在60-80歲、射血分?jǐn)?shù)為35%、血肌酐正常的患者，用上述列線圖進(jìn)行預(yù)測的總評(píng)分為118分（其中年齡54.5分，射血分?jǐn)?shù)為36分，血肌酐為27.5分）。該患者（118分）存活到100天的概率大約為81%，存活半年概率為70%；具有這些特征的患者，生存率為75%時(shí)的生存時(shí)間為至少140天（大約）。

注：Q25（25%分位數(shù)）生存時(shí)間為什么對(duì)應(yīng)的是生存率是75%而不是25%？這里是否可以說“具有這些特征的患者，75%的人能存活的時(shí)間至少是140天”？對(duì)此有疑問的可以參見此次筆記最后的【困惑】部分。

##下面看下單獨(dú)預(yù)測分位數(shù)函數(shù)的fun.side設(shè)置

Nom.Quan.psm<-nomogram(HFDeath.psm,lp=F,fun=list(Q25,median), fun.at=c(seq(20,100,by=20),seq(150,300,by=50),seq(400,1200,by=200),1800,2500), funlabel=c("Q25 Survival Time","Q50 Survival Time"))

plot(Nom.Quan.psm,fun.side=c(rep(1,12),3,1,1,1),xfrac=.25,col.grid = gray(c(0.8, 0.95)))

參考

rms包幫助文件：https://cran./web/packages/rms/

Zhongheng Zhang,Michael W. Kattan.Drawing Nomograms with R: applications to categorical outcome and survival data.Ann Transl Med 2017;5(10):211. [http://dx./10.21037/atm.2017.04.01]

【困惑的鑰匙】

最后說一下預(yù)測分位數(shù)生存時(shí)間帶給我的困惑。

上圖中的“Q50生存時(shí)間”在rms幫助文件中是用中位生存時(shí)間（MST，Median Survival Time）來描述的。從定義上來講，分位數(shù)是針對(duì)一個(gè)數(shù)據(jù)集而言的，比如25分位數(shù)Q25，表示在這個(gè)數(shù)據(jù)集中有25%的數(shù)據(jù)小于Q25這數(shù)。但在生存分析中，第Q分位數(shù)生存時(shí)間（Qth Quantile Survival Time）并不是“所有生存時(shí)間的Q分位數(shù)”。本例“所有生存時(shí)間的25%、50%、75%分位數(shù)”為73天、115天、205天。一個(gè)生存料在隨訪結(jié)束時(shí)，所有受試者的生存時(shí)間可以計(jì)算獲得一個(gè)中位值（本例115天），但是可能因?yàn)槔鄯e生存率尚未降低到50%而不會(huì)存在中位生存時(shí)間（本例就未達(dá)到中位生存時(shí)間），根本原因就是累積生存率的計(jì)算考慮到了刪失，刪失的存在可能會(huì)讓Q分位數(shù)的生存時(shí)間比“所有生存時(shí)間的Q分位數(shù)”要長的多。而且筆者感覺對(duì)這個(gè)分位數(shù)生存時(shí)間的定義在不同的軟件中可能不太一樣。比如在R的survival包中，明確說“生存曲線S (t)的第k分位數(shù)是一條高度為p=1-k的水平線與S (t)曲線圖相交的位置[The kth quantile for a survival curve S(t) is the location at which a horizontal line at height p=1-k intersects the plot of S(t)]”。對(duì)應(yīng)的，第Q百分位數(shù)生存時(shí)間是描述的是S(t)≤1-Q%時(shí)所對(duì)應(yīng)的生存時(shí)間，在這種定義下你會(huì)發(fā)現(xiàn)，越小的百分位生存時(shí)間其對(duì)應(yīng)的生存時(shí)間越短。但在SPSS中卻是百分位數(shù)越小其對(duì)應(yīng)的生存時(shí)間越長，顯然這里對(duì)第Q百分位數(shù)生存時(shí)間的定義是生存率S(t)為Q%時(shí)生存曲線對(duì)應(yīng)的生存時(shí)間。當(dāng)然中位生存時(shí)間處于中間的位置，不論哪種定義結(jié)果都是一樣的。另外，由于刪失的存在，一般采用乘積極限法來估算累積生存率，第Q分位數(shù)的生存時(shí)間和(1-Q%)的患者尚存活的時(shí)間可能不一致，所以嚴(yán)格地講，中位生存時(shí)間是指累積生存率為50%時(shí)所對(duì)應(yīng)的時(shí)間，而不是存活人數(shù)尚有一半（死亡人數(shù)達(dá)到一半）時(shí)的時(shí)間。

不論怎么講，分位數(shù)生存時(shí)間是針對(duì)一個(gè)觀察樣本而不是單獨(dú)某一個(gè)個(gè)體而言的。這樣就有一個(gè)疑問，當(dāng)我們用列線圖來預(yù)測時(shí)，預(yù)測對(duì)象是某一個(gè)受試者，一個(gè)受試者怎么會(huì)有的分位數(shù)生存時(shí)間呢？而且分位數(shù)不應(yīng)該是一個(gè)值嗎，為什么列線圖中是一個(gè)范圍？

rms包的幫助文件中對(duì)Quantile()的描述為“The Quantile method for cph returns an S function for computing quantiles of survival time (median, by default)”，“Quantile.psm create S functions that evaluate the quantile functions analytically, as functions of time and the linear predictor values”。我甚至一度認(rèn)為幫助文件中出現(xiàn)的“中位生存時(shí)間（Median Survival Time）”是“生存時(shí)間的中位值”。也曾認(rèn)為應(yīng)該像分位數(shù)回歸那樣來理解這個(gè)范圍，把預(yù)測軸前面的分位數(shù)生存時(shí)間（因變量）作為解讀的一個(gè)條件，后面的軸刻度是預(yù)測的患者生存時(shí)間，如“生存率75%時(shí)仍然存活患者中，具有這些特征的患者可以存活到140天”。

實(shí)際上不僅這個(gè)分位數(shù)生存時(shí)間，在某個(gè)隨訪時(shí)間點(diǎn)，生存概率、平均生存時(shí)間其實(shí)都應(yīng)該是一個(gè)值，我們?cè)陬A(yù)測軸上看到的也是一個(gè)范圍。這其實(shí)并不難理解，因?yàn)橐粋€(gè)特征組合（解釋變量）是一個(gè)值，不同的特征組合就對(duì)應(yīng)著不同的預(yù)測值了。至于為什么一個(gè)人會(huì)有平均生存時(shí)間、分位數(shù)生存時(shí)間，用理解生存概率的辦法來理解這幾個(gè)概念會(huì)更容易一些。但只有一個(gè)人是沒法計(jì)算概率的，對(duì)一個(gè)人而言，要么生要么死，不可能既生又死（比如某人生存概率60%就有點(diǎn)這個(gè)既生又死的味道）。其實(shí)當(dāng)預(yù)測模型告訴你“小A的1年生存概率是60%”的時(shí)候，指的是擁有像小A特征的那些人1年生存概率是60%，而小A剛好是其中一個(gè)而已。同樣的，對(duì)于分位數(shù)生存時(shí)間、平均生存時(shí)間也是如此，真正的含義應(yīng)該是具有某些特征組合的這一類人，其某個(gè)分位數(shù)的生存時(shí)間和平均生存時(shí)間。

本例隨訪結(jié)束，累積生存率為57.6%，尚未到達(dá)中位生存時(shí)間，因此在SPSS里面是不會(huì)結(jié)算這個(gè)值，如果預(yù)測模型使用中位生存時(shí)間和75分位數(shù)生存時(shí)間屬于“超綱”預(yù)測。對(duì)于這種數(shù)據(jù)，可以考慮使用25分位數(shù)生存時(shí)間（Q25ST=100）（SPSS中是75分位數(shù)生存時(shí)間）來描述。