上一篇关注单因素多因素cox模型的构建并在lung数据集中进行了实战
survival包学习笔记——cox回归（一） - 简书 (jianshu.com)
接下来我们继续上一篇谈一谈如何评价我们构建的模型，并且进一步可视化回归结果。
一、函数介绍
上篇我们已经介绍了cox回归的一大函数，survival包中的coxph函数，今天我们来说说另一自产自销的函数，rms::cph
cph(formula = formula(data), data=environment(formula),
weights, subset, na.action=na.delete,
method=c("efron","breslow","exact","model.frame","model.matrix"),
singular.ok=FALSE, robust=FALSE,
model=FALSE, x=FALSE, y=FALSE, se.fit=FALSE,
linear.predictors=TRUE, residuals=TRUE, nonames=FALSE,
eps=1e-4, init, iter.max=10, tol=1e-9, surv=FALSE, time.inc,
type=NULL, vartype=NULL, debug=FALSE, ...)
formula: 公式与coxph一致，使用Surv包装一下。Surv(time, status) ~ 变量
data：数据集
x：默认值为 FALSE。设置为 TRUE 可将扩展的design matrix作为返回的拟合对象的元素 x
y：默认值为 FALSE。设置为 TRUE 以返回响应值的向量（Surv 对象）作为拟合的元素 y
二、如何评价模型？
1. 区分度评价指标：C指数（C-Index），重新分类指数（Net reclassification index，NRI）
2. 一致性评价指标：校正曲线（Calibration plot）
3. 临床有效性评价指标：决策分析曲线（Decision Curve Analysis， DCA）
C指数: survival::coxph函数构建模型后自动计算，校正曲线和DCA都是通过rms构建的
1、Concordance index：
在使用survival::coxph时，C指数倒是不需要再求，summary一下拟合模型就可以看到：
或者
一般而言，我们将0.51-0.7认为是低可信度，0.71-0.9为中等可信度，> 0.9为高可信度
2、校正曲线（Calibration plot）
【参考】
临床预测模型列线图的评估：校准曲线 - 简书 (jianshu.com)
calibration_curve(校准曲线): 分类模型可视化技术之一 - 知乎 (zhihu.com)
模型构建完成后需要对模型进行评估和验证其性能。模型预测的生存率与实际的差距有多大呢？一般是看校准曲线。
例：一个模型（其C指数为0.8）评估某位患者5年复发率为70%。说明该模型有80%的把握确认复发率=70%。那70%这个数与实际相差有多大呢，那就需要看校准曲线了。
从这个例子可以看出，C指数或AUC值是判断模型的区分能力的，即有多大把握预测复发率为70%，而校准曲线是看与预测与实际相符程度的，即预测的这个70%复发率与实际复发率有多大差别。
3、决策分析曲线（Decision Curve Analysis， DCA）
【参考】
决策曲线分析法（Decision Curve Analysis，DCA）曲线 | Public Library of Bioinformatics (plob.org)
手把手教你用R画列线图(Nomogram)及解读结果 - 知乎 (zhihu.com)
三、cox回归结果展示
Cox结果主要包括：HR, HR_95%low, HR_95%high, P
我们最常用的展示方式就是森林图，同时森林图也是我们展示cox回归的最佳方式，我们来看发表的森林图：
森林图绘制
数据准备，将cox回归的结果读入，
rm(list = ls())
load("RData/7. Clinical/MSC cox result.RData")
#森林图unicox_OS----
df <- unicox_OS[unicox_OS$p_val<0.1,]
rownames(df)[1] <- "MSC1"
gene <- rownames(df)
colnames(df)
hr <- sprintf("%.3f",df$"HR")
hrLow <- sprintf("%.3f",df$"HR_l")
hrHigh <- sprintf("%.3f",df$"HR_H")
Hazard.ratio <- paste0(hr," (",hrLow,"-",hrHigh,")")
pVal <- ifelse(df$p_val<0.01, "<0.01", sprintf("%.3f", df$p_val))
使用基础函数绘制森林图，大致思路：将画布按照3：2分为两部分，将字和线段逐个打印上
#森林图
if(T){
pdf(file="unicox_OS.pdf", width = 8,height = 3.5)
n <- nrow(df)
nRow <- n+1
ylim <- c(1-0.5,nRow+0.5)
layout(matrix(c(1,2),nc=2),width=c(3,2))
#绘制森林图左边的基因信息
par(mar=c(4,0.5,1,0))#上左下右
xlim = c(0.4,2.5)
plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,xlab="",ylab="")
text.cex=1
text(0.4,n:1,gene,adj=0,cex=text.cex)
text(1.3,n:1,pVal,adj=1,cex=text.cex)
text(2.5,n:1,Hazard.ratio,adj=1,cex=text.cex)
text(1.35,n+1,'P-value',cex=text.cex,font=2,adj=1)
text(2.35,n+1,'HR (95% CI)',cex=text.cex,font=2,adj=1,)
abline(h=n+1.3,col="black",lty=1,lwd=2);
abline(h=n+0.6,col="black",lty=1,lwd=2);
abline(h=0.6,col="black",lty=1,lwd=2)
#绘制森林图
par(mar=c(4,0,2,1),mgp=c(2,0.5,0))
xlim = c(0,3)
plot(1,xlim=xlim,ylim=ylim,type="n",axes=F,ylab="",xaxs="i",xlab="")
arrows(as.numeric(hrLow),n:1,
as.numeric(hrHigh),n:1,
angle=0,code=3,length=0.01,col="black",lwd=2.5)
abline(v=1.0,col="black",lty=2,lwd=1.5)
boxcolor ="black"
points(as.numeric(hr), n:1, pch = 18, cex=2,col = boxcolor)
axis(1)
dev.off()
}
当然，还有不少包装好的R包可以直接绘制森林图，感兴趣的小伙伴可以探索一下
总之呐，图是为了更好的展示，因此，结果最重要，对吧
