上一篇關(guān)注單因素多因素cox模型的構(gòu)建并在lung數(shù)據(jù)集中進行了實戰(zhàn)
survival包學(xué)習(xí)筆記——cox回歸（一） - 簡書 (jianshu.com)
接下來我們繼續(xù)上一篇談一談如何評價我們構(gòu)建的模型，并且進一步可視化回歸結(jié)果蚂夕。
一来吩、函數(shù)介紹
上篇我們已經(jīng)介紹了cox回歸的一大函數(shù)凉倚，survival包中的coxph函數(shù)，今天我們來說說另一自產(chǎn)自銷的函數(shù)买置，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：數(shù)據(jù)集
x：默認值為 FALSE。設(shè)置為 TRUE 可將擴展的design matrix作為返回的擬合對象的元素 x
y：默認值為 FALSE蒸苇。設(shè)置為 TRUE 以返回響應(yīng)值的向量（Surv 對象）作為擬合的元素 y
二磷蛹、如何評價模型？
1. 區(qū)分度評價指標：C指數(shù)（C-Index）溪烤，重新分類指數(shù)（Net reclassification index味咳，NRI）
2. 一致性評價指標：校正曲線（Calibration plot）
3. 臨床有效性評價指標：決策分析曲線（Decision Curve Analysis庇勃， DCA）
C指數(shù): survival::coxph函數(shù)構(gòu)建模型后自動計算，校正曲線和DCA都是通過rms構(gòu)建的
1槽驶、Concordance index：
在使用survival::coxph時责嚷，C指數(shù)倒是不需要再求，summary一下擬合模型就可以看到：
或者
一般而言掂铐，我們將0.51-0.7認為是低可信度罕拂，0.71-0.9為中等可信度，> 0.9為高可信度
2全陨、校正曲線（Calibration plot）
【參考】
臨床預(yù)測模型列線圖的評估：校準曲線 - 簡書 (jianshu.com)
calibration_curve(校準曲線): 分類模型可視化技術(shù)之一 - 知乎 (zhihu.com)
模型構(gòu)建完成后需要對模型進行評估和驗證其性能爆班。模型預(yù)測的生存率與實際的差距有多大呢？一般是看校準曲線辱姨。
例：一個模型（其C指數(shù)為0.8）評估某位患者5年復(fù)發(fā)率為70%柿菩。說明該模型有80%的把握確認復(fù)發(fā)率=70%。那70%這個數(shù)與實際相差有多大呢炮叶，那就需要看校準曲線了碗旅。
從這個例子可以看出，C指數(shù)或AUC值是判斷模型的區(qū)分能力的镜悉，即有多大把握預(yù)測復(fù)發(fā)率為70%祟辟，而校準曲線是看與預(yù)測與實際相符程度的，即預(yù)測的這個70%復(fù)發(fā)率與實際復(fù)發(fā)率有多大差別侣肄。
3旧困、決策分析曲線（Decision Curve Analysis， DCA）
【參考】
決策曲線分析法（Decision Curve Analysis稼锅，DCA）曲線 | Public Library of Bioinformatics (plob.org)
手把手教你用R畫列線圖(Nomogram)及解讀結(jié)果 - 知乎 (zhihu.com)
三吼具、cox回歸結(jié)果展示
Cox結(jié)果主要包括：HR, HR_95%low, HR_95%high, P
我們最常用的展示方式就是森林圖，同時森林圖也是我們展示cox回歸的最佳方式矩距，我們來看發(fā)表的森林圖：
森林圖繪制
數(shù)據(jù)準備拗盒，將cox回歸的結(jié)果讀入，
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))
使用基礎(chǔ)函數(shù)繪制森林圖锥债，大致思路：將畫布按照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()
}
當(dāng)然，還有不少包裝好的R包可以直接繪制森林圖哮肚，感興趣的小伙伴可以探索一下
總之吶登夫，圖是為了更好的展示，因此允趟，結(jié)果最重要恼策，對吧
