建立一个模型后,我们常常会评价模型的区分度(discrimination)和校准度(calibration)。生存模型中,我们经常会看到使用calibration图来呈现模型的校准度。笔者近期查阅了网上许多绘制calibration图的R代码,发现很多代码忽略了time.inc参数的使用,甚至注释中提到了这个参数的重要性却仍然用错了。将愚见记录与此,若有理解错误,烦请大佬指正。
发现这个问题是源于丁香园中一个讨论(关于R语言做calibration curve的相关问题),其中一个人提出“需要注意的是cph()建立模型的时候,time.inc参数的数值应该和calibrate()函数中参数u保持一致,例如3年就都是36。”
笔者首先查阅网上的相关代码,有的代码中会注释到time.inc需与参数u一致,但是却将time.inc局限地理解为建立Nomogram时的利用predictor计算survival的生存函数中的time参数。实际上在建立cph()时就应该预先设定好相应的time.inc参数的值。之前笔者也不确定这个参数的设置是否真的这么重要,然后就在自己的数据中进行尝试和比较,发现得到的calibration图差别比较大,于是找到rms包的源代码,发现其中calibrate.cph.s文件中如下几行和这个问题相关:
#获取生存预测模型的summary
ssum <- fit$surv.summary
if(!length(ssum)) stop('did not use surv=TRUE for cph( )')
cat("Using Cox survival estimates at ", dimnames(ssum)[[1]][2],
" ", unit, "s\n", sep="")
#获取每个层的2nd time(即u时刻)的平均生存率S #注意:此处不是基线风险函数h0
surv.by.strata <- ssum[2, , 1] #2nd time= at u, all strata
xb <- fit$linear.predictors #获取Cox模型的predictor #注意:这里线性部分是减去了均值的
if(length(stra <- fit$strata))
surv.by.strata <- surv.by.strata[stra]
#根据公式计算每个样本的生存率
survival <- as.vector(surv.by.strata ^ exp(xb))
其中我们看到,计算的基线生存函数,是u时刻的基线生存函数,但是这个u值却直接取自fit的模型的summary中,并非我们之后在调用calibrate方法时设置的u参数(这么做也是为了使整体源代码更简洁,但却容易造成用户误用)。
根据这些,笔者认为,正确的calibration应该这么绘制(以下举例,且忽略了nomogram的绘制)
coxm <- cph(Surv(time, status==1)~Age.+Grade.+Stage.+Surgery.+Size.+Site.,
x=T,y=T,data=All.data, surv=T, time.inc=36)
cal<-calibrate(coxm,u=36,cmethod='KM',m=150, B=100)
plot(cal,xlim = c(0,1),ylim= c(0,1),
errbar.col=c(rgb(0,0,0,maxColorValue=255)),
col=c(rgb(255,0,0,maxColorValue=255)),
xlab = "Predicated 3-year CSS",
ylab = "Actual 3-year CSS")
abline(0,1,lty=3,lwd=2,col=c(rgb(0,0,255,maxColorValue= 255)))
以上愚见,欢迎大家一起进一步讨论这个问题!