【原】临床预测模型[序]：来看看这样建立的模型哪里需要改进？

前几天在PubMed检索了一下“Prediction Model”，然后就看到了预测模型在每年的发表情况，真的可以用指数增长来描述其爆发式增长了。

精准医学已经是被谈了多年的一个概念了，精准医学所追求的极致应该是个体化医疗，如果能算命般的准确预测某个患者的临床结局或者某种干预对结局的影响，这的确是一件了不起的事情。临床预测模型（Clinical Prediction Model）就是利用已知特征通过模型来预测未知。临床预测模型一般建立在多因素回归分析（如线性回归、logistic回归、Cox回归等）的基础上，是用预测变量（自变量）来进行疾病诊断、疾病是否发生或者疾病的预后情况的预测，据此可分为诊断模型、预后模型及疾病预测模型。

打算花一段时间从头到尾地演示一下临床预测模型的建立与评估。

在开始建立一个模型建立前，你最好先明确一下你的目的：你建立的这个模型是用来干啥的（临床问题）？在设计好你的临床方案（试验设计、预测变量、临床结局、样本量、统计模型）后，开始进行数据收集。获得数据后我们可能还需要花很长的一段时间进行数据清理，并进行数据的预处理（缺失值处理、预测因素异常值处理、预测因素编码(连续变量？分类变量？样条？)）。其后是预测变量的筛选（专业确定？先单后多？逐步回归？最优子集？Lasso回归？），回归模型建立，并以一定呈现形式来展示模型的结果。模型建立以后还有一个重要的步骤，就是对模型进行内部评估和外部评估（曲ROC线下面积(AUC、C指数)、校准曲线）。如果建立了多个模型，则还会涉及模型间的比较（净重分类指数NRI、综合区分改善指数IDI、临床决策曲线DCA等）。

这是一个大工程，需要分多次笔记来完成。

示例：SMART研究（Second Manifestations of ARTerial disease study）。该研究是一项由荷兰乌得勒支大学医学中心协调的正在进行的前瞻性队列研究。此次示例的数据纳入了1996年9月至2006年3月期间纳入研究的3873名患者，纳入的患者均有动脉粥样硬化的临床表现（短暂性脑缺血发作、缺血性卒中、外周动脉疾病、腹主动脉瘤或冠心病），知情同意后进入研究。研究终点事件是（急性）血管性死亡、（非）致死性缺血性卒中或（非）致死性心肌梗死，以及任何这些血管事件的复合终点。在心血管疾病领域的许多预测模型都是根据没有动脉粥样硬化临床表现的受试者数据开发的，本研究则是建立一个针对已存在心血管疾病（冠状动脉疾病、脑动脉疾病、外周动脉疾病和腹主动脉瘤）患者的预测模型，用于估计血管事件（中风、心肌梗死或心血管死亡）1年、3年和5年的发生风险。

示例来源：Ewout W.Steyerberg.Clinical Prediction Models: A Practical Approach to Development, Validation, and Updating(Second Edition,2019).

涉及的潜在解释变量有人口学特征、经典的危险因素、血压水平、血脂水平、既往动脉粥样硬化史及动脉粥样硬化标志物。具体有：stdays（生存时间，天），event（终点事件，0=删失，1=发生）；Sex（性别，1=male，2=female），Age（年龄，岁）；BMI（体质指数）,Smoking（吸烟，0=Never,1=Former,2=Current），Packyears（包年数,(每天吸烟量/20)*吸烟年数），Alcohol（饮酒，0=Never,1=Former,2=Current），Diabetes（糖尿病，0=No,1=Yes）；SBP（收缩压,mmHg），DBP（舒张压,mmHg）；TC（总胆固醇,mmol/L），HDLC（高密度脂蛋白胆固醇,mmol/L），LDLC（低密度脂蛋白胆固醇,mmol/L），TG（甘油三酯,mmol/L）；Cerebral（脑血管病史，0=No,1=Yes），Coronary（心血管病史，0=No,1=Yes），Peripheral（外周血管病史，0=No,1=Yes），AAA（腹主动脉瘤病史，0=No,1=Yes）；Hcy（同型半胱氨酸,µmol/L），Gln（谷氨酰胺,µmol/L），CrCl（肌酐清除率,mL/min），AU（白蛋白尿，0=No,1=Micro,2=Macro）,IMT（内膜中层厚度,mm）、CAStn（颈动脉狭窄>50%，0=No,1=Yes）。

你可能觉得这简单呀，不就是建立一个Cox生存模型然后画个列线图，最后评估一下就可以了吗？我们不妨先这样做一下看。采用的软件是R4.1.3。

#导入数据，并将数据集设置为数据框

library(haven)

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

smart<-smartdata[c("stdays","event","Sex","Age","BMI","Smoking","Alcohol","Diabetes","MAP","LDLC","HDLC","TG","Cerebral","Coronary","Peripheral","AAA","Hcy","Gln","CrCl","AU","IMT","CAStn")]##为什么选择要分析的变量重新建立一个数据库，而不是直接通过命令smart<-as.data.frame(as.matrix(smardatat))来实现达到目的？主要是考虑到有些没有被选择的变量可能存在缺失值，直接通过前面这个命令可能会损失更多的信息

smart<-as.data.frame(as.matrix(smart))

#将分类变量设置为因子。二分类变量其实可以不用设置分类变量，但为了后面绘制列线图中更具有可读性（不是0和1而是0和1代表的具体意义），可将分类变量全部设置为因子

smart$Sex<-factor(smart$Sex,levels=c(1,2),labels=c("male","female"))

smart$Smoking<-factor(smart$Smoking,levels=c(0,1,2),labels=c("Never","Former","Current"))

smart$Alcohol<-factor(smart$Alcohol,levels=c(0,1,2),labels=c("Never","Former","Current"))

smart$Diabetes<-factor(smart$Diabetes,levels=c(0,1),labels=c("No","Yes"))

smart$Cerebral<-factor(smart$Cerebral,levels=c(0,1),labels=c("No","Yes"))

smart$Coronary<-factor(smart$Coronary,levels=c(0,1),labels=c("No","Yes"))

smart$Peripheral<-factor(smart$Peripheral,levels=c(0,1),labels=c("No","Yes"))

smart$AAA<-factor(smart$AAA,levels=c(0,1),labels=c("No","Yes"))

smart$AU<-factor(smart$AU,levels=c(0,1,2),labels=c("No","Micro","Macro"))

smart$CAStn<-factor(smart$CAStn,levels=c(0,1),labels=c("No","Yes"))

#删除缺失值记录。大部分函数对无法分析有缺失的记录，甚至报错

smart<-na.omit(smart) 

 示例共有病例3873例，其中460例发生失效事件，相对而言预测变量并不多，因此考虑不进行先单后多的初步筛选，而是想直接采用全子集回归，但是20个变量的全子集回归的计算量还是超出了我的预估：电脑跑了一晚上竟然没跑出来，最终采用了向后逐步回归，筛选标准为AIC。考虑到血压中SBP与DBP存在较强的相关性，可以只选取对心血管预测有更重要作用的SBP，或者是同时考虑了SBP和DBP的平均动脉压MAP（DBP+(SBP-DBP)/3）。另外，血脂水平中TC与LDLC等存在较强的相关性，目前专业上也普遍认为TC或LDLC是心脑血管疾病独立的危险因素，而HDLC则可能是保护因素，因此模型对血脂方面的预测因素考虑LDLC、HDLC和TG。

#逐步回归（向后法）

library(StepReg)

stepwiseCox(Surv(stdays,event==1)~Sex+Age+BMI+Smoking+Alcohol+Diabetes+MAP+LDLC+HDLC+TG+Cerebral+Coronary+Peripheral+AAA+Hcy+Gln+CrCl+AU+IMT+CAStn,data=smart, selection ="backward", select ="AIC",method=c("efron"))

如上图结果所示，以AIC为标准，向后逐步回归筛选到10个变量（如果P为标准筛选到的是9个变量）。

下面用函数nomogram{rms}以这10个变量绘制列线图。

library(rms)

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

CoxNM<-datadist(smart)

options(datadist='CoxNM')

Coxfit<-cph(Surv(stdays,event==1)~Age+Alcohol+Cerebral+Coronary+Peripheral+AAA+Hcy+Gln+CrCl+IMT,data=smart,x=T,y=T,surv=T)

surv<-Survival(Coxfit)  ##获得预测生存概率

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

surv3<-function(x)surv(3*365.25,lp=x)

surv5<-function(x)surv(5*365.25,lp=x)

NomCox<-nomogram(Coxfit,fun=list(surv1,surv3,surv5),lp=F,funlabel=c("1-Year Survival","3-Year Survival","5-Year Survival"),maxscale=100,fun.at=c('0.95','0.85','0.70','0.6','0.5','0.4','0.3','0.2','0.1','0.05','0.01'))

plot((NomCox))

你也可以利用regplot函数将列线图具有交互性，看起来更漂亮一些。我们在《绘制更有颜值的列线图》中有介绍。

library(survival)

Coxfit<-coxph(Surv(stdays,event==1)~Age+Alcohol+Cerebral+Coronary+Peripheral+AAA+Hcy+Gln+CrCl+IMT,data=smart)

library(regplot)

regplot(Coxfit, plots=c("density","boxes"), observation=T, title="Survival Nomogram", failtime=c(365.25,1095.75,1826.25), prfail=F, clickable=TRUE, points=TRUE, dencol="green",  boxcol="yellow", droplines=TRUE)

模型建立完毕，我们还需要对模型验证一下：结果显示C-index=0.708，表明模型区分度一般般；校准曲线与y=x的直线重合性很差，表明模型的校准度也不好，即结局的实际发生概率与模型的预测概率一致性并不好。

#区分度C-Statistic

CI_Coxfit<-coxph(Surv(stdays,event==1)~Age+Alcohol+Cerebral+Coronary+Peripheral+AAA+Hcy+Gln+CrCl+IMT,x=T,y=T,data=smart)

sum.cox<-summary(CI_Coxfit)

CI<-sum.cox$concordance

CI

         C      se(C)

0.70843993 0.02530677

#校准曲线Calibration Curve

CalCox<-calibrate(Coxfit,cmethod='KM',method='boot',u=5*365.25,m=300,B=200)

plot(CalCox,xlim=c(0,1.0),ylim=c(0,1.0))

为了更清楚的显示，我们放大一下局部：

plot(CalCox,xlim=c(0.9,1.0),ylim=c(0.5,1.0))

估计会有不少人也会像上面这样做。

这样做有什么问题吗？

<可以思考一下的>

……

……

……

你可能还会用不同的方法来筛选变量，比如向前逐步回归、双向逐步回归、lasso回归，然后建立不同的模型并进行比较选择比较好的一个。这个确实值得去做，既然没有上帝之眼，就多做几个选个最好的吧。

我们常常忽略的还有缺失值和异常值的处理。

本例实际上存在着大量的缺失值。总病例虽然有3873例（460例失效事件），但如果只分析完整记录的病例，在去掉潜在的20个变量有缺失值的记录后，仅剩1953例（其中失效事件133例，删失1820例），根据10倍EPV（Events Per Variable，每个变量对应10个阳性事件。当结果阳性事件更多时则是阴性事件满足10EPV）的经验原则，似乎能分析的变量也不是想象中的那么多，更为关键的是直接的暴力删除丢失了大量的信息。怎么办？这就需要进行缺失值的插补。

异常值也是一个非常值得关注的重要问题，因为异常值可能会对回归模型造成重大影响。将连续型变量调整为分类变量是一种选择。在本示例出处的教材(Clinical Prediction Models(Second Edition,2019))中，Ewout W Steyerberg给出了简单易操作的方法，值得借鉴：Biologically implausible values were set to missing values, and remaining extreme values were winsorized by shifting the values approximately below the 1st centile and above the 99th centile to"truncation points".简单地说，对生物学上不可能出现的值进行删除，其他的一些极值采用缩尾值进行处理，所谓缩尾处理就是将超出变量特定百分位范围的数值替换为特定百分位数值的方法，而不是删掉。文中给出的是将低于第1百分位和高于第99百分位为“截断点”。

还有在建模的时候，变量以什么样的形式纳入到模型中更合适？是连续变量还是分类变量形式？需不需要进行变换？需不需要将具有相似效果的预测因素进行合并？本例在分析的时候仅简单涉及但其实差得还很远。

本例最后对模型的评估用的是内部验证，如果我们还有一个验证集，外部验证也是很多人想知道的……

……

总之问题还有很多，实在很难用一篇笔记就能面面俱到，此次笔记就算作一个引子吧，计划后面再分专题记录这其中的林林总总。

— END —