(二)临床模型评价:
当临床预测模型成功建立,就会面临一个问题,如何去评价模型的质量。模型的评价一般是通过比较人群的模型预测结果与实际观测结果,来评价预测模型的效果。
对于模型的评价指标,我们往往分成三种:
1. 区分度评价指标:C指数(C-Index),重新分类指数(Net reclassification index,NRI);
2. 一致性评价指标:校正曲线(Calibration plot);
3. 临床有效性评价指标:决策分析曲线(Decision Curve Analysis, DCA);
今天就让阿琛来给大家分别讲讲如何评价前期建立的临床预测模型。
1.建立预测模型:
首先,我们一起快速的回顾一下如何构建Cox回归模型。
1.1 引用R包:
#install.packages("rms")
#install.packages("foreign")
#install.packages("survival")
library(rms)
library(foreign)
library(survival)
1.2 读取文件:
setwd("C:Users 00Desktop13_clinical") #设置工作目录
rt <- read.table("clinical.txt",header=T,sep=" ") #读取数据
head(rt) #查看前5行的数据
在该数据集中,共包含6个不同的临床病理特征以及患者的生存状态和生存时间。
1.3 构建Cox回归模型:
首先,对数据的相关参数进行整理与设置:
rt$gender <- factor(rt$gender,labels=c("F", "M"))
rt$stage <- factor(rt$stage,labels=c("Stage1", "Stage2", "Stage3", "Stage4"))
rt$T <- factor(rt$T,labels=c("T1", "T2", "T3", "T4"))
rt$M <- factor(rt$M,labels=c("M0", "M1"))
rt$N <- factor(rt$N,labels=c("N0", "N1", "N2", "N3"))
rt$risk <- factor(rt$risk,labels=c("low", "high"))
str(rt)
ddist <- datadist(rt) #将数据打包
options(datadist='ddist')
units(rt$futime) <- "year" #定义时间的单位
接着,将所有变量纳入模型中,构建Cox回归模型:
#构建回归模型
f <- cph(Surv(futime, fustat) ~gender + stage +T + M + N + risk,
x=T, y=T, surv=T, data=rt)
f #查看模型的相关参数
自此,我们的Cox回归模型就构建完成了。接下来就是从不同的角度来评价模型的质量。
2.C指数(C-index)的计算:
validate(f, method="boot", B=1000, dxy=T) #重复模拟1000次
rcorrcens(Surv(futime, fustat) ~ predict(f), data = rt)
结果显示:
在此模型中,C-index为0.664,即1-0.336=0.664。
当然,根据C指数而言,该模型的预测区分度并不是太高。一般而言,我们将0.51-0.7认为是低可信度,0.71-0.9为中等可信度,> 0.9为高可信度。
3.Calibration plot的绘制:
绘制三年生存率的校准曲线
#计算三年生存率的校准曲线
f3 <- cph(Surv(futime, fustat) ~gender + stage +T + M + N + risk,
x=T, y=T, surv=T, time.inc = 3, data=rt)
cal3 <- calibrate(f3, cmethod="KM", method="boot", u= 3, m= 90, B= 1000)
#一般将患者分成3组,即m约等于患者总数/3,且重复模拟1000次
plot(cal3)
结果显示:
绘制五年生存率的校准曲线
#计算五年生存率的校准曲线
f5 <- cph(Surv(futime, fustat) ~gender + stage +T + M + N + risk,
x=T, y=T, surv=T, time.inc = 5, data=rt)
cal5 <- calibrate(f5, cmethod="KM", method="boot", u= 5, m= 90, B= 1000)
plot(cal5)
结果显示:
对于Calibration plot,横坐标为模型预测得到的患者生存率,纵坐标为患者实际观察得到的生存率,三点之间的联系与虚线之间越相近,越能说明模型良好的一致性。
4.DCA曲线的绘制:
source("stdca.R") #把stdca.R放在之前设定的工作目录中,并进行引用
rt$three.years.Survival.Probabilitynew = c(1- (summary(survfit(f, newdata = rt),times = 3)$surv)) #计算3年生存率
rt$five.years.Survival.Probabilitynew = c(1- (summary(survfit(f, newdata = rt),times = 5)$surv)) #计算5年生存率
write.csv(rt, "DCA.Survival.csv")
#将患者的3年和5年生存率结果保存
随后,根据患者的生存情况,分别绘制3年和5年的DCA曲线;
#绘制3年生存率的DCA曲线
stdca(data=rt,
outcome="fustat", ttoutcome="futime", timepoint=3,
predictors="three.years.Survival.Probabilitynew",
xstop=0.9, smooth=TRUE)
结果显示:
#绘制5年生存率的DCA曲线
stdca(data=rt,
outcome="fustat", ttoutcome="futime", timepoint=5,
predictors="five.years.Survival.Probabilitynew",
xstop=0.9,smooth=TRUE)
结果显示:
5.NRI指数的计算:
对于NRI指数的计算,有2个专门的R包,分别为nricens和PredictABEL;其中,nricens计算出来的是绝对的NRI,PredictABEL则为相对的NRI,一般在使用过程中选择其中之一即可。
#install.packages("nricens")
library(nricens)
随后,分别构建两个不同的模型,以进行比较:
mstd <- cph(Surv(futime, fustat) ~gender + stage +T + M + N,
x=T, y=T, surv=T, data=rt) #使用除risk外的变量进行建模
mnew <- cph(Surv(futime, fustat) ~gender + stage +T + M + N + risk,
x=T, y=T, surv=T, data=rt) #该模型包含所有的变量
接着,计算NRI指数;
#计算连续的NRI,取值为5%时就定义为重分类
#3年生存率NRI
nricens(mdl.std = mstd, mdl.new = mnew, t0 = 3, updown = 'diff',cut = 0.05, niter = 200)
#5年生存率NRI
nricens(mdl.std = mstd, mdl.new = mnew, t0 = 5, updown = 'diff',cut = 0.05, niter = 200)
#计算分类变量计算的NRI
nricens(mdl.std = mstd, mdl.new = mnew, t0 = 5, updown = 'category',cut = c(0.3,0.6), niter = 200)
对于NRI的计算,其主要是比较新旧两个模型之间存在区分度。其中,updown主要定义了一个样本的风险是否变动的方式,category是指分类值,即熟悉的低、中、高风险,另有一种diff,为连续值。
Cut是判断风险高低的临界值。当updown为diff时,cut只需设置1个值,比如0.05,即认为当预测的风险在新旧模型中相差5%时,即被认为是重新分类了;而当updown为category时,可以对cut设置两个不同的值,即0~29%为低风险,30%~59%为中风险,60%~100%为高风险。
结果显示:
当updown为diff时,3年和5年的NRI值分别为0.26和0.53,均大于5%,可以认为预测的风险在新旧模型中发生了重新分类;
3年生存率NRI的95%可信区间:
5年生存率NRI的95%可信区间:
同时,由于做了200次重复取样,相当于有200个NRI,于是NRI就有了标准误和95%的置信区间。同理,可以得到当updown为category时的NRI值。
到此,对于临床预测模型评价指标的介绍就到此结束了。大家可以对照着整个分析过程,使用示例数据或者自己的模型数据,进行相应的学习~