R|可视化|multiROC

作者：理查德·高

1.数据准备

• futime生存时间：删除time<30d、unknown和空值；

• fustat生存状态：删除unknown、空值；

• age：删除unknown、空值；

• sex：删除unknown、空值，female=0，male=1；

• stage：stage I=1，stage II=2，stage III =3, stage IV =4;

• T分期：T1a/T1b/...=1, T2a/T2b/...=2, T3a/T3b/...=3, T4a/T4b/...=4

• N分期：N0=0, N1=1,N2=2,N3=3,去掉Nx和unknown

• M分期：M0=0，M1=1,去掉Mx和unknown

• grade：grade1=1，grade2=2，grade3=3，grade4=4，去掉Gx和unknown

indepInput

2.单因素预后

rm(list=ls())
library(survival)
library(survminer)
rt <- read.table("indepInput.txt",header=T,sep="\t",check.names=F,row.names=1)
outTab=data.frame()
for(i in colnames(rt[,3:ncol(rt)])){
    x=coxph(Surv(futime,fustat)~rt[,i],data=rt)
    x=summary(x)
    outTab=rbind(outTab,
                cbind(id=i,
                     HR=signif(x$conf.int[,"exp(coef)"],2),
                     HR.95L=signif(x$conf.int[,"lower .95"],3),
                     HR.95H=signif(x$conf.int[,"upper .95"],3),
                     pvalue=x$coefficients[,"Pr(>|z|)"]))
}
outTab[,-1] <- apply(outTab[,-1],2,as.numeric)
outTab$fdr <- p.adjust(outTab$pvalue,method="fdr",length(outTab$pvalue))
#write.table(outTab,file="uniCox.txt",sep="\t",row.names=F,quote=F)

select1 <- c("futime","fustat",outTab$id[outTab$pvalue<0.05])#提取显著的character
indepInput1 <- indepInput[,select1]
save(indepInput1,file="multicoxInput.Rdata")

3.多因素预后

rm(list=ls())
library(survival)
library(survminer)
load("multicoxInput.Rdata")
multiCox=coxph(Surv(futime,fustat)~.,data=indepInput1)
x = summary(multiCox)
outTab = data.frame()
outTab = cbind(coeff=x$coefficients[,1],
              HR=x$conf.int[,"exp(coef)"],
              HR.95L=x$conf.int[,"lower .95"],
              HR.95H=x$conf.int[,"upper .95"],
              pvalue=x$coefficients[,"Pr(>|z|)"])
outTab=cbind(id=rownames(outTab),outTab)
#write.table(outTab,file="multiCox.txt",sep="\t",row.names=F,quote=F)

4. multiROC

rm(list = ls())
library(survivalROC)
rt <- read.table("./indepInput.txt",sep = "\t",row.names = 1,header = T,check.names = F)
pdf(file = "./multiROC.pdf",width = 6,height = 6)
rt$futime <- rt$futime/365
rocCol <- rainbow(ncol(rt)-2)
aucText=c()

par(oma=c(0.5,1,0,1),font.lab=1.5,font.axis=1.5)
roc=survivalROC(Stime = rt$futime,status = rt$fustat,marker = rt$riskScore,predict.time = 1,method = "KM")
plot(roc$FP,roc$TP,
     type = "l",xlim = c(0,1),ylim = c(0,1),col=rocCol[1],
     xlab="1-Specificity", ylab="Sensitivity",
     lwd=2,cex.main=1.3,cex.lab=1.2,cex.axis=1.2,font=1.2)
aucText=c(aucText,paste0("risk score"," (AUC=",sprintf("%.3f",roc$AUC),")"))
abline(0,1)

j=1
for (i in colnames(rt[,3:(ncol(rt)-1)])){
  roc=survivalROC(Stime = rt$futime,status = rt$fustat,marker = rt[,i],predict.time = 1,method = "KM")
  j=j+1
  aucText=c(aucText,paste0(i," (AUC=",sprintf("%.3f",roc$AUC),")"))
  lines(roc$FP,roc$TP,type = "l",xlim = c(0,1),ylim = c(0,1),col=rocCol[j],lwd=2)
}
legend("bottomright",aucText,lwd = 2,bty = "n",col = rocCol)
dev.off()

multiROC

参考文献：

survivalROC使用手册