R语言潜类别分析 r语言cindex

roc曲线和cindex曲线属于一致性的指标。后台有粉丝问怎么绘制cox回归中的时间相关性roc曲线和时间相关性cindex曲线，今天来演示一下。

我们先导入数据和R包

library(survival)
library("survminer")
library(foreign)
bc <- read.spss("E:/r/test/Breast cancer survival agec.sav",
                use.value.labels=F, to.data.frame=T)
bc <- na.omit(bc)
names(bc)

这是我们常用的乳腺癌数据，（公众号回复：乳腺癌,可以获得该数据），age表示年龄，pathsize表示病理肿瘤大小（厘米），lnpos表示腋窝淋巴结阳性，histgrad表示病理组织学等级，er表示雌激素受体状态，pr表示孕激素受体状态，status结局事件是否死亡，pathscat表示病理肿瘤大小类别（分组变量），ln_yesno表示是否有淋巴结肿大，time是生存时间，后面的agec是我们自己设定的，不用管它。

有部分变量为分类变量，我们先把它转换成因子

bc$histgrad<-as.factor(bc$histgrad)
bc$er<-as.factor(bc$er)
bc$pr<-as.factor(bc$pr)
bc$ln_yesno<-as.factor(bc$ln_yesno)

我们先建立3个cox归回方程，随便建的

f1<-coxph(Surv(time,status)~er+histgrad+pr+age+ln_yesno,bc,x=TRUE,y=TRUE)
f2<-coxph(Surv(time,status)~er+histgrad+ln_yesno,bc,x=TRUE,y=TRUE)
f3<-coxph(Surv(time,status)~ln_yesno,bc,x=TRUE,y=TRUE)

很多R包可以做时间相关性ROC,本次我们介绍riskRegression包，如大家感兴趣下次可以介绍一下其他的包

library(riskRegression)

假设我们想了解f1的时间相关性roc, 这里要注意一下formula=Surv(time,status)~1等于是不含权重的，如果把1改成协变量等于进行加权

A3<- riskRegression::Score(list("f1"=f1),
                           formula=Surv(time,status)~1,
                           data=bc,
                           metrics="auc",
                           null.model=F,
                           times=seq(3,132,1))
plotAUC(A3)

这样图形就生成了，非常简单，我们还可以把数据提取出来放到ggplot上绘图

auc<-plotAUC(A3)
ggplot()+geom_line(data=auc, aes(times,AUC),linetype=1,size=1,alpha = 0.6,colour="red")+
  geom_ribbon(data=auc, aes(times,ymin = lower, ymax = upper),alpha = 0.1,fill="red")+
  geom_hline(yintercept=1, linetype=2,size=1)+theme_classic()+ 
  labs(title = "时间相关性ROC", x="times", y="AUC")

假设我们想了解f1和f2的时间相关性roc

A3<- riskRegression::Score(list("f1"=f1,"f2"=f2),
                           formula=Surv(time,status)~1,
                           data=bc,
                           metrics="AUC",
                           null.model=F,
                           times=seq(3,132,1))
auc<-plotAUC(A3)
ggplot()+geom_line(data=auc, aes(times,AUC,group=model,col=model))+
  geom_ribbon(data=auc, aes(times,ymin = lower, ymax = upper,fill=model),alpha = 0.1)+
  geom_hline(yintercept=1, linetype=2,size=1)+theme_classic()+ 
  labs(title = "时间相关性ROC", x="times", y="AUC")

下面来介绍一下时间先关性cindex，这个需要使用到pec包，先画f1的

library(pec)
A1<-pec::cindex(list("f1"=f1),
                formula=Surv(time,status)~er+histgrad+pr+age+ln_yesno,
                data=bc,
                eval.times=seq(3,132,1))
plot(A1)

假设想了解f1、f、f3先关性cindex

A1<-pec::cindex(("f1"=f1,"f2"=f2,"f3"=f3),
                formula=Surv(time,status)~er+histgrad+pr+age+ln_yesno,
                data=bc,
                eval.times=seq(3,132,1))
plot(A1)

还可以在颜色，标签等进行更改，我这里就不弄了。