转自个人微信公粽号【易学统计】的统计学习笔记：R语言实现Cox模型校准度曲线绘制
研究背景
这是关于cox模型的第二篇文章，上一篇文章分享了运用Lasso回归如何筛选变量，将筛选后的变量绘制Nomogram图，本章分享构建模型后，如何绘制校准曲线。
cox模型的验证不同于Logistic回归，cox的结局包括时间和状态，所以对于某个患者来说，他的结果是否准确，就要看模型在他随访的时间点，所预测的结局是否和真实的状态一致，两者一致说明拟合不错，模型准确。虽然大部分患者的患者的随访时间都不一样，但对于生存分析来说，我们更关注一段时间在群体中的效应，比如1年生存率、3年生存率和5年生存率，以及中位生存率。换句话说，在评价手术后或者药物治疗后的效果时，更关注某一段时间后的平均疗效。
案例研究
本文数据采用R自带的lung数据集，一个晚期肺癌数据集，收集了228例癌症患者的生存资料，包含患者年龄、性别、生存时间、生存状态等10个变量。本文利用年龄、性别和ecog评分三个变量构建模型，绘制列线图并绘制指定时间点的校准曲线。
临床研究一般有提供多个危险因素，首先做单因素的筛选，具体筛选方法，见公众号之前的文章。详细Nomogram图绘制见之前文章，本章详细说明校准曲线的绘制和参数说明，采用验证的方法是基于原始数据集的重抽样验证。
R代码及解读
##加载包 明确每个包的作用
library(rms) ##绘制列线图
library(survival) ##生存分析包
##调用数据，数据格式与普通的spss中格式一样，一行代表一条观测，
##一列代表一个变量；
data(lung)
d <- lung
str(d)
aggr(d,prop=T,numbers=T) #判断数据缺失情况，红色表示有缺失。
d <- na.omit(d) #按行删除缺失值
str(dt) ##查看每个变量结构
aggr(dt,prop=T,numbers=T) #判断数据缺失情况，红色表示有缺失。
d <- na.omit(dt) 按行删除缺失值
由图片可看到所有变量都为蓝色，没有缺失值，红色表示有缺失，并展示出缺失比例。如果用na.omit()函数按照行删除。
第一步，数据整理。
###添加变量标签，在列线图上展示分类标签，作图用到
d$sex <- factor(d$sex,levels = c(1,2),labels = c('male','female'))
##########################################################################
d[,6] <- as.factor(d[,6])
##结局变量转换
d$status <- ifelse(d$status==2,1,0)
str(d)##可以查看数据结构
第二步：构建模型并绘制列线图
dd<-datadist(d) #设置工作环境变量，将数据整合
options(datadist='dd') #设置工作环境变量，将数据整合
##
coxm <- cph(Surv(time,status)~age+sex+ph.ecog,x=T,y=T,data=d,surv=T)
###绘制cox回归生存概率的nomogram图
## 构建Nomo图的对象只能是rms保重d额cph()函数
nom <- nomogram(coxm,fun=list(function(x)surv(1*365,x), ##算出不同时间节点生存率值
function(x)surv(2*365,x)), ##算出不同时间节点生存率值，显示在列线图上
funlabel = c('1-year probability',
'2-year probability'),
lp=F,
fun.at=c('0.9','0.85','0.80','0.70','0.6','0.5','0.4','0.3','0.2','0.1'))
par(mar=c(2,5,3,2),cex=0.8)##mar 图形空白边界 cex 文本和符号大小
plot(nom,xfrac=0.6)
该图为1年和2年生存率的列线图，性别的标签在图上展示出来了。接下来分别对1年生存率和2年生存率做验证曲线。。
第三步：绘制校准曲线
###这里演示，采用age sex 和ph.ecog 构建coxph模型
##构建校准曲线
##time.in 和 u 要是一样的，都是要评价的时间节点
coxm_1 <- cph(Surv(time,status)~age+sex+ph.ecog,data=d,surv=T,x=T,y=T,time.inc = 365)
cal_1<-calibrate(coxm_1,u=365,cmethod='KM',m=30,B=1000)
##绘制1年生存期校准曲线
par(mar=c(7,4,4,3),cex=1.0)
plot(cal_1,lwd=2,lty=1, ##设置线条形状和尺寸
errbar.col=c(rgb(0,118,192,maxColorValue = 255)), ##设置一个颜色
xlab='Nomogram-Predicted Probability of 1-year DFS',#便签
ylab='Actual 1-year DFS(proportion)',#标签
col=c(rgb(192,98,83,maxColorValue = 255)),#设置一个颜色
xlim = c(0,1),ylim = c(0,1)) ##x轴和y轴范围
##绘制2年生存期校曲线
##time.in 和 u 要是一样的，都是要评价的时间节点
coxm_2 <- cph(Surv(time,status)~age+sex+ph.ecog,data=d,surv=T,x=T,y=T,time.inc = 2*365)
cal_2<-calibrate(coxm_2,u=2*365,cmethod='KM',m=30,B=1000)
plot(cal_2,lwd=2,lty=1, ##设置线条宽度和线条类型
errbar.col=c(rgb(0,118,192,maxColorValue = 255)), ##设置一个颜色
xlab='Nomogram-Predicted Probability of 1-year DFS',#便签
ylab='Actual 2-year DFS(proportion)',#标签
col=c(rgb(192,98,83,maxColorValue = 255)),#设置一个颜色
xlim = c(0,1),ylim = c(0,1)) ##x轴和y轴范围
绘制校准曲线的参数说明：
1.绘制校准曲线前需要在模型函数中添加参数x=T，y=T和time.inc 参数
2.u需要和之前模型中定义好的time.inc一致，根据原数据，要以天为单位
即365或者365/2，u和time.inc其实是一个意思，表示要评价的时间节点。如果它俩不一样，表示你想要评价的时间点不一样，矫正的结果当然不同了。time.inc全称是 time increase 即时间增量。当surv=T时候，表示要计算surv .summary数组，此时time.inc用于计算数组内部不同时间点的数据(seq(0,maxtime,by=time.inc))。评价时间点不能小于最大时间的1/25，也不能超过最大随访时间，会超出评价评价限度,不能评价。如果time.inc省略，则会根据时间单位自动选择，这往往不是我们要的，所以还需要指定值。
3.m要根据样本量来确定，由于校准曲线一般将所有的样本均分为5组(对应图中的5分节点)，而m代表每组样本量数，因此m*5等于或者近似等于样本量
4.b代表最大再抽样的样本量，一般b=1000，足够使用。
5.由两个校准图可看，横坐标为预测的生存率，纵坐标为实际的生存率，对角线是预测概率等于实际概率，偏离对角线越远说明预测的误差越大。以上两个图中预测曲线基本和对角线重合，结果还行，这一点也可以通过计算C-index指数看出端倪。
以上就是Cox模型校准曲线的所有内容了。后面还有更多高分统计方法分享，请持续关注哦~
如果您觉得有用，请点赞，转发哦~
更多统计小知识，请关看 公粽号 易学统计
更多阅读
R语言|中位生存时间列线图绘制
基于Lasso回归筛选变量构建Cox模型并绘制Nomogram
R语言Logistic回归模型验证及Nomogram绘制
如何进行高维变量筛选和特征选择(一)？Lasso回归
