生存资料（参数、半参数）预测模型的列线图绘制

列线图（Alignment Diagram），也称诺莫图(Nomogram），绝是临床预测模型最为常见的一种呈现形式。列线图可以根据受试者的一些特征（解释变量）来预测结局事件的发生概率或者患者的生存率等。生存资料预测模型的列线图除了可以预测生存率，还可以预测分位数生存时间、生存时间期望均值、风险比。

生存资料预测模型的构建绝大多数是基于COX回归法。COX回归属于半参数法，在多因素的生存分析中几乎是霸主般的存在，主要是因为COX回归几乎不需要考虑前提条件，仅要求满足比例风险假定就可以了。除了COX回归，参数生存模型也是生存资料非常重要的分析方法。但参数法需要数据满足一定分布，一说到分布啥的一些理论，基本上就阻断了大部分人继续下去的欲望，甚至头皮发麻恶心想吐，哪还有心思去搞清楚怎么用。但其实如果满足适用条件，参数法可能会比非参数法更好地解释数据。正如在比较组间效应时，如果数据满足独立正态同方差，我们更愿意适用t检验/方差分析，而不会首选非参数检验一样，虽然使用非参数检验也没有错。与生存时间相关的一些重要分布有"weibull", "exponential", "gaussian", "logistic","lognormal" and "loglogistic"。

此次笔记的任务就是绘制生存资料预测模型的列线图，包括COX回归模型和参数生存模型。

数据来源：UCI Machine Learning Repository

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

这是一个心力衰竭临床记录数据集，有13个临床特征。笔者对个别字段的数据类型进行了转换。

下面直接上命令清单：

导入数据与因子设置

##导入数据

library(haven)

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

HF<-as.data.frame(HF)

##设置因子及因子各水平的标签

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回归模型

在经过单变量分析和全子集回归的分析后（此处略），最终选择了三个变量（Age.cat、EF、SCr.cat）来建立COX生存模型，模型正态基本满足PH假定。

1.2.1 COX回归模型列线图绘制

##载入列线图绘制程序包

library(rms)

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

##设置工作环境，汇总预测变量分布信息。该设置为必需项。

CoxNM<-datadist(HF)

options(datadist='CoxNM')

##设置预测变量在列线图上显示的名称。不设置则默认为变量的名称

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：结（ties）的处理方式有efron、breslow和exact，默认为efron。这里改为method="breslow"是为了得到跟SPSS一致的结果，STATA默认的也是breslow。同样的数据在采用R的coxph{survival}和cph{rms}的进行Cox回归，其结果与SPSS有略微差异，如果你像我一样纠结，可以在R中把结的处理方式改一下就可以了；

#x、y、surv：这几个参数是为了绘制后面的列线图方便，都设置为TRUE就可以了；

#time.inc：这在后期进行模型评估中可能会用到，比如评估模型在不同时间点上的校准度，常常用函数calibrate绘制1年、3、5年的校准曲线，这时候就需要用time.inc指定时间点，然后再calibrate函数的u参数进行呼应。

##建立生存率的估计函数

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函数可以对cph创建的对象返回一个S函数，用于计算生存率的估计值。除了Survival函数，Quantile()返回的S函数可用于计算生存时间的分位数（默认为中位数，可通过q来指定具体的分位数）；Mean()返回返回一个计算平均生存时间的函数。参数估计模型还有Hazard()用于估计风险函数；

#通过function函数定义了三个生存率的计算函数，这些函数以一个时间点和线性预测器作为参数，并返回该时间点的生存率。本例时间单位是天，用91.31天和182.62天来表示1个月和3个月，1年应该用365.25天，但本例最长的随访时间是285天，如果使用大于285的数字，结果会报错“Error in approx(fu[s], xseq[s], fat, ties = mean) :调用内插至少需要两个非NA值的数”。此处仅示例。

##列线图的参数设置

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."))

##列线图绘制的主要函数，参数又多又复杂又难理解，仅介绍几个常用的、重要的

#lp：默认设置为FALSE，即不显示线性预测器（linear predictor）的那条坐标轴。我们最终想预测的概率就是这个线性预测器的值通过一定的函数变换得来的。列线图的基本原理就是根据预测变量对结局的影响大小（回归系数），给各个预测变量的水平进行赋分，评分之和通过函数转化成结局事件发生概率。这个线性预测值可以理解为承接变量总评分与预测概率之间的中间变量，对列线图的使用而言其实没啥毛用；

#fun：对线性预测值进行转换的函数；

#fun.at：预测值横轴的刻度，本例是在坐标轴上显示的概率值。可使用该参数的默认值，也可通过该参数改成需要显示的值。同一个fun函数可以是使用同一个fun.at向量值，不同fun函数需要用fun.at值的向量列表标识。本例有三个自建函数，不过都是对应生成生存率的Survival这一个函数，因此三条预测概率坐标轴刻度使用了相同的fun.at值;

#funlabel：预测横轴名称，默认"Predicted Value"；

#conf.int: 是否显示得分的置信区间，默认否；

#maxscale：默认单变量最大分数值是100分。

##绘制、美化列线图

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

##主要是将前面Nomogram函数的结果画出来，一些参数可以起到美化作用

#lplabel：线性预测器横轴的名称，默认"Linear Predictor"。不过这个变量我一般会在Nomogram中关掉不显示；

#fun.side：预测轴上的刻度值位置，1用于定位在轴下方（默认值），3用于轴上方。预测横轴的刻度值有时候会非常拥挤甚至重叠在一起，如一个轴有5个标记标签值，第二个和第三个相互重叠，可该参数进行指定：fun.side=c(1,1,3,1,1)；

#col.conf：置信区间的颜色设置；

#xfrac：各刻度线横轴与轴名称之间的空白间隔所占的比例，默认0.35。如果你觉得坐标轴名称与坐标轴隔得距离太远，可以把这个值设置的小一些；

#cex.axis：刻度轴刻度标签的字符大小，默认0.85；

#cex.var：刻度轴名称的字符大小，默认1；刻度轴标题默认是变量名称，可通过label()进行修改；

#col.grid ：默认不绘制垂直参考网格线。该参数若长度为1，则对主要和次要参考线使用相同的颜色，长度为2则分别对应主要和次要参考线的颜色，常用灰色col.grid = gray(c(0.8, 0.95))；

#points.label：单变量评分轴的名称，默认为'Points'；

#total.points.label：总评分轴的名称，默认'Total Points'。

列线图的主要控制参数及绘制结果如下，列线图的解读方法可参见第2部分参数生存模型的列线图部分。

1.2.2带有交互作用的列线图绘制

模型是否存在交互作用应在1.2步进行考察，此处仅演示其列线图的绘制操作。

##上接设置工作环境变量显示名称

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)

#除了建立生存率预测函数，还建立一个生存时间期望均值的预测函数。最后一次随访时结局为非删失事件，Mean函数会有报警。由于删失的存在，平均生存时间是一个有偏估计量，一般用生存曲线下面积或累计发病率曲线上面积表示。

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参数的：fun.side=list(rep(1,10),rep(3,10),rep(c(1,3),5),c(rep(1,10),3,1,3,1))，结果报错。单独预测生存率或者单独预测生存时间均值时都没有问题，但合在一起就出现了错误，寻因良久，未果，只得作罢。不过可以适当对fun.at参数进行调整来解决刻度值拥挤的问题。

对某一个患者，在预测时交互作用只选择其中符合该患者特征那一行进行评分即可，比如要预测的患者70岁，EF和Age.cat的交互作用我们只用列线图EF(Age.cat=60-80)这一个横轴，结合EF值来打分就可以了。列线图的使用可参见第2部分参数生存模型的列线图部分。

【2】参数生存模型

rms包中的psm函数可以提供6种参数生存模型（Parametric Survival Model）的拟合，包括威布尔分布（weibull），指数分布（exponential）、正态分布（gaussian）、logistic、对数正态（lognormal）和对数logistic。

经初步考察，在psm函数中提供的6个参数分布中，数据拟合对数正态分布似乎更好一些。

##上接设置工作环境变量显示名称

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)

#除了建立生存率预测函数，同时建立生存时间分位数的预测函数。默认是中位值，通过q参数可指定任意分位数

#本例25分位数生存时间为100天，随访结束时的累积生存率为57.6%，尚未到达中位生存时间

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"))

#采用默认的fun.at，分位数预测刻度轴的标签会非常得拥挤，但在调整展示两个函数（此处是生存函数和分为数函数）刻度值的位置时，fun.side参数再次报错，只能使用fun.at可以进行适当的调整。

一个年龄在60-80岁、射血分数为35%、血肌酐正常的患者，用上述列线图进行预测的总评分为118分（其中年龄54.5分，射血分数为36分，血肌酐为27.5分）。该患者（118分）存活到100天的概率大约为81%，存活半年概率为70%；具有这些特征的患者，生存率为75%时的生存时间为至少140天（大约）。

注：Q25（25%分位数）生存时间为什么对应的是生存率是75%而不是25%？这里是否可以说“具有这些特征的患者，75%的人能存活的时间至少是140天”？对此有疑问的可以参见此次笔记最后的【困惑】部分。

##下面看下单独预测分位数函数的fun.side设置

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.r-project.org/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.doi.org/10.21037/atm.2017.04.01]

【困惑的钥匙】

最后说一下预测分位数生存时间带给我的困惑。

上图中的“Q50生存时间”在rms帮助文件中是用中位生存时间（MST，Median Survival Time）来描述的。从定义上来讲，分位数是针对一个数据集而言的，比如25分位数Q25，表示在这个数据集中有25%的数据小于Q25这数。但在生存分析中，第Q分位数生存时间（Qth Quantile Survival Time）并不是“所有生存时间的Q分位数”。本例“所有生存时间的25%、50%、75%分位数”为73天、115天、205天。一个生存料在随访结束时，所有受试者的生存时间可以计算获得一个中位值（本例115天），但是可能因为累积生存率尚未降低到50%而不会存在中位生存时间（本例就未达到中位生存时间），根本原因就是累积生存率的计算考虑到了删失，删失的存在可能会让Q分位数的生存时间比“所有生存时间的Q分位数”要长的多。而且笔者感觉对这个分位数生存时间的定义在不同的软件中可能不太一样。比如在R的survival包中，明确说“生存曲线S (t)的第k分位数是一条高度为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)]”。对应的，第Q百分位数生存时间是描述的是S(t)≤1-Q%时所对应的生存时间，在这种定义下你会发现，越小的百分位生存时间其对应的生存时间越短。但在SPSS中却是百分位数越小其对应的生存时间越长，显然这里对第Q百分位数生存时间的定义是生存率S(t)为Q%时生存曲线对应的生存时间。当然中位生存时间处于中间的位置，不论哪种定义结果都是一样的。另外，由于删失的存在，一般采用乘积极限法来估算累积生存率，第Q分位数的生存时间和(1-Q%)的患者尚存活的时间可能不一致，所以严格地讲，中位生存时间是指累积生存率为50%时所对应的时间，而不是存活人数尚有一半（死亡人数达到一半）时的时间。

不论怎么讲，分位数生存时间是针对一个观察样本而不是单独某一个个体而言的。这样就有一个疑问，当我们用列线图来预测时，预测对象是某一个受试者，一个受试者怎么会有的分位数生存时间呢？而且分位数不应该是一个值吗，为什么列线图中是一个范围？

rms包的帮助文件中对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”。我甚至一度认为帮助文件中出现的“中位生存时间（Median Survival Time）”是“生存时间的中位值”。也曾认为应该像分位数回归那样来理解这个范围，把预测轴前面的分位数生存时间（因变量）作为解读的一个条件，后面的轴刻度是预测的患者生存时间，如“生存率75%时仍然存活患者中，具有这些特征的患者可以存活到140天”。

实际上不仅这个分位数生存时间，在某个随访时间点，生存概率、平均生存时间其实都应该是一个值，我们在预测轴上看到的也是一个范围。这其实并不难理解，因为一个特征组合（解释变量）是一个值，不同的特征组合就对应着不同的预测值了。至于为什么一个人会有平均生存时间、分位数生存时间，用理解生存概率的办法来理解这几个概念会更容易一些。但只有一个人是没法计算概率的，对一个人而言，要么生要么死，不可能既生又死（比如某人生存概率60%就有点这个既生又死的味道）。其实当预测模型告诉你“小A的1年生存概率是60%”的时候，指的是拥有像小A特征的那些人1年生存概率是60%，而小A刚好是其中一个而已。同样的，对于分位数生存时间、平均生存时间也是如此，真正的含义应该是具有某些特征组合的这一类人，其某个分位数的生存时间和平均生存时间。

本例随访结束，累积生存率为57.6%，尚未到达中位生存时间，因此在SPSS里面是不会结算这个值，如果预测模型使用中位生存时间和75分位数生存时间属于“超纲”预测。对于这种数据，可以考虑使用25分位数生存时间（Q25ST=100）（SPSS中是75分位数生存时间）来描述。

END