分享一段比较不同模型的C-index和RMS的代码

已封装成函数,如下: compareRiskModal=function(os,ev,risks){ library('survRM2') library('Hmisc') require(survival) library(survcomp) library(prodlim) os=data$SURVIVAL_...

已封装成函数,如下:

compareRiskModal=function(os,ev,risks,max_time){
  library('survRM2')
  library('Hmisc')
  require(survival)
  library(survcomp)
  library(prodlim)
  
  # Parameters for RMS
  tau <-max_time # Restriction point
  n.grid <-100 # For plotting purposes
  # Set up data using age as the marker
  OS <-Surv(os,ev)
  marker.pp<-seq(from=0,to=1,length=n.grid)
  mdat=cbind()
  c.all=list()
  c.indexs=c()
  for(i in 1:ncol(risks)){
    marker <-as.numeric(risks[,i])
    marker.qq<-quantile(marker,marker.pp)
    fitdat.df<-data.frame(marker=marker)
    newdat.df<-data.frame(marker=marker.qq)
    # Calculations
    cox.model1<-coxph(OS~marker,data=fitdat.df)
    rms.calc <-summary(survfit(cox.model1,newdata=newdat.df),rmean=tau)
    rms.mean <-rms.calc$table[,"*rmean"]
    mdat=cbind(mdat,rms.mean)
    c1=concordance.index(x=as.numeric(risks[,i]), surv.time=os, surv.event=ev,
                      method="noether")
    c.all=c(c.all,list(c1))
    c.indexs=c(c.indexs,paste0(colnames(risks)[i],':',round(c1$c.index,2),',95%CI(',round(c1$lower,2),'-',round(c1$upper,2),'),p=',signif(c1$p.value,2)))
  }
  
  mt=ceiling(max(mdat))
  # RMS Curve
  cls=rainbow(ncol(risks))
  if(ncol(risks)>1){
    plot(marker.pp,as.numeric(mdat[,1]),pch=20,col = cls[1], cex=0.5,ylim=c(0,mt),xlab="",ylab="")
    if(ncol(risks)>2){
      for(i in 2:(ncol(risks)-1)){
        par(new=TRUE)
        plot(marker.pp,as.numeric(mdat[,i]),pch=20,col = cls[i], cex=0.5,ylim=c(0,mt),xlab="",ylab="")  
      }
    }
    par(new=TRUE)
    plot(marker.pp,as.numeric(mdat[,ncol(risks)]),pch=20,col = cls[ncol(risks)]
         , cex=0.5,ylim=c(0,mt),xlab="Percentile of Score",ylab="RMS")
  }else{
    plot(marker.pp,as.numeric(mdat[,1]),pch=20,col = "red", cex=0.5,ylim=c(0,mt),xlab="Percentile of Score",ylab="RMS")
  }
  legend("topright",   c.indexs, pch=20,col=cls, title= "C-index", inset = .05,cex = 0.7)
  #print(c.all)
  if(ncol(risks)>1){
    all.cp=c()
  for(i in 1:(ncol(risks)-1)){
    c1=concordance.index(x=risks[,i], surv.time=os, surv.event=ev,
                         method="noether")
    for(j in (i+1):ncol(risks)){
      c2=concordance.index(x=risks[,j], surv.time=os, surv.event=ev,
                           method="noether")
      p=min(cindex.comp(c1, c2)$p.value,cindex.comp(c2, c1)$p.value)
      all.cp=c(all.cp,paste0(colnames(risks)[i],'-vs-',colnames(risks)[j],':p=',signif(p,2)))
    }
    legend("bottomleft",   all.cp, title= "Compare C-index", inset = .05,cex = 0.7)
  }
  }
}

compareRiskModal(data$SURVIVAL_TIM,data$STATUS,data.frame(IPV=data$IPVHCC,MOD=data$Risk),5*365)

调用规则:

os:生存时间

ev: 生存事件

risks:一个数据框,列为各个模型的风险得分,行为样本和os,ev需一一对应。

max_time:计算的最大时间,比如五年 

体验效果如:

attachments-2019-11-9SGHwRdt5dca5c7c6cbe7.png


  • 发表于 2019-11-12 15:16
  • 阅读 ( 4330 )
  • 分类:编程语言

相关问题

1 条评论

请先 登录 后评论
不写代码的码农
祝让飞

生物信息工程师

118 篇文章

作家榜 »

  1. 祝让飞 118 文章
  2. 柚子 91 文章
  3. 刘永鑫 64 文章
  4. admin 57 文章
  5. 生信分析流 55 文章
  6. SXR 44 文章
  7. 张海伦 31 文章
  8. 爽儿 25 文章