一、繪制符合ggplot2風格的圖片，可以加theme

1、先定義一個函數(shù)，生成timeROC對象，注意數(shù)據(jù)集和相應列名需要修改

library(survivalROC)
## Define a helper functio nto evaluate at various t
survivalROC_helper <- function(t) {
    survivalROC(Stime        = ovarian$futime,
                status       = ovarian$fustat,
                marker       = ovarian$lp,
                predict.time = t,
                method       = "NNE",
                span = 0.25 * nrow(ovarian)^(-0.20))
}

2、計算每180天的ROC參數(shù)，具體時間可以修改，傳入數(shù)據(jù)為向量

## Evaluate every 180 days
survivalROC_data <- data_frame(t = 180 * c(1,2,3,4,5,6)) %>%
    mutate(survivalROC = map(t, survivalROC_helper),
           ## Extract scalar AUC
           auc = map_dbl(survivalROC, magrittr::extract2, "AUC"),
           ## Put cut off dependent values in a data_frame
           df_survivalROC = map(survivalROC, function(obj) {
               as_data_frame(obj[c("cut.values","TP","FP")])
           })) %>%
    dplyr::select(-survivalROC) %>%
    unnest() %>%
    arrange(t, FP, TP)

3、畫圖

## Plot
survivalROC_data %>%
    ggplot(mapping = aes(x = FP, y = TP)) +
    geom_point() +
    geom_line() +
    geom_label(data = survivalROC_data %>% dplyr::select(t,auc) %>% unique,
               mapping = aes(label = sprintf("%.3f", auc)), x = 0.5, y = 0.5) +
    facet_wrap( ~ t) +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5),
          legend.key = element_blank(),
          plot.title = element_text(hjust = 0.5),
          strip.background = element_blank())

二、用timeROC包畫，便于事后各種對比

library(timeROC)
time_roc_1 <-  timeROC(
  T = data_ti_roc$STIATime,
  delta = data_ti_roc$STIA,
  marker = data_ti_roc$lp_ti_roc_1, #方向相反加個-
  cause = 1,
  weighting="marginal", #uses the Kaplan-Meier
  times = seq(10,60,10),
  ROC = TRUE,
  iid = T
)

library(RColorBrewer)  # 使用包之前，先加載
plot(time_roc_1,time=60,title =FALSE,col = brewer.pal(4,'Dark2')[1])        
plot(time_roc_2,time=60,add=TRUE,col = brewer.pal(4,'Dark2')[2]) 
plot(time_roc_3,time=60,add=TRUE,col = brewer.pal(4,'Dark2')[3]) 
plot(time_roc_4,time=60,add=TRUE,col = brewer.pal(4,'Dark2')[4])
legend("bottomright",c("CHA2DS2VASC","CHA2DS2VASC+PTFV15000","CHA2DS2VASC+LAE","CHA2DS2VASC+PTFV15000+LAE"),col = brewer.pal(4,'Dark2'),lty=1,lwd=2)

三、繪制timeAUC比較曲線

帶置信區(qū)間的并不好看

## 繪制timeAUC
par(oma = c(0,0,0,0),      #外框據(jù)整個繪圖區(qū)域四周距離，單位是"行"
    mar = c(4,4,0.5,0.5),  #圖片四周距外框距離,同上，如果用omi或mri,單位是英寸，不算四周文字
    mgp = c(1.7,0.5,0),    #三個數(shù)字，分別為坐標標題到坐標軸距離，坐標數(shù)字到軸距離,ticks到軸距離
    cex = 0.8            #字體大小，不包括lengend
    #tck = -0.01           #ticks大小及方向，負值向外
)
plotAUCcurve(time_roc_1,conf.int=F,col="#00468B99")
legend(x = 20,y=1,c("CHA2DS2VASC"),col=c("#00468B99"),lty=1,lwd=2,cex = 0.8)
plotAUCcurve(time_roc_4,conf.int=F,col="#ED000099",add = T)
legend(x = 20,y=0.95,c("CHA2DS2VASC+PTFV15000+LAE"),col=c("#ED000099"),lty=1,lwd=2,cex = 0.8)

四、繪制timeAUC比較圖

##------------------------------------------繪制timeAUC比較---------------------------------------
pdf(file = "time_auc.pdf",width = 4.0,height = 4.0)
## 繪制timeAUC
par(oma = c(0,0,0,0),      #外框據(jù)整個繪圖區(qū)域四周距離，單位是"行"
    mar = c(4,4,0.5,0.5),  #圖片四周距外框距離,同上，如果用omi或mri,單位是英寸，不算四周文字
    mgp = c(1.7,0.5,0),    #三個數(shù)字，分別為坐標標題到坐標軸距離，坐標數(shù)字到軸距離,ticks到軸距離
    cex = 0.8            #字體大小，不包括lengend
    #tck = -0.01           #ticks大小及方向，負值向外
)
plotAUCcurve(time_roc_1,conf.int=F,col="#00468B99")
legend(x = 20,y=1,c("CHA2DS2VASC"),col=c("#00468B99"),lty=1,lwd=2,cex = 0.8)
plotAUCcurve(time_roc_4,conf.int=F,col="#ED000099",add = T)
legend(x = 20,y=0.95,c("CHA2DS2VASC+PTFV15000+LAE"),col=c("#ED000099"),lty=1,lwd=2,cex = 0.8)
dev.off()

五、timeROC面積比較

##---------------------------------timeROC各項參比較----------------------------------------

ci_roc_1 = confint(time_roc_1)
ci_roc_2 = confint(time_roc_2)
ci_roc_3 = confint(time_roc_3)
ci_roc_4 = confint(time_roc_4)

ci_roc_1 <- ci_roc_1$CI_AUC %>% t() %>% data.frame() %>% dplyr::select(t.60)
ci_roc_2 <- ci_roc_2$CI_AUC %>% t() %>% data.frame() %>% dplyr::select(t.60)
ci_roc_3 <- ci_roc_3$CI_AUC %>% t() %>% data.frame() %>% dplyr::select(t.60)
ci_roc_4 <- ci_roc_4$CI_AUC %>% t() %>% data.frame() %>% dplyr::select(t.60)

df_ti_rocs <- data.frame(AUC = c(time_roc_1$AUC['t=60'],time_roc_2$AUC['t=60'],time_roc_3$AUC['t=60'],time_roc_4$AUC['t=60']),
                         lower = c(ci_roc_1[1,1],ci_roc_2[1,1],ci_roc_3[1,1],ci_roc_4[1,1]),
                         upper = c(ci_roc_1[2,1],ci_roc_2[2,1],ci_roc_3[2,1],ci_roc_4[2,1])
)

df_ti_rocs$ORCI <- df_ti_rocs %>% apply(1,function(x){str_c(c(round(x["AUC"],3),"(",round(x["lower"],4),"-",round(x["upper"],4),")"),collapse = "")})

df_ti_rocs <- df_ti_rocs %>% mutate(p = c(NA,
                                          compare(time_roc_2,time_roc_1)$p_values_AUC[[6]],
                                          compare(time_roc_3,time_roc_1)$p_values_AUC[[6]],
                                          compare(time_roc_4,time_roc_1)$p_values_AUC[[6]]
))

write.csv(df_ti_rocs,file = "auc_compare.csv")

六、計算cox模型的Harrell's C-Statistics

##----------------------------------------------計算cox模型的Harrell's  C-Statistics------------------------------------
#BiocManager::install("survcomp")
library(survcomp)

cindex_with_ci <- fit_ti_rocs %>% lapply(function(x){c_index <- concordance.index(predict(x,data),
                                                                                  surv.time = data$STIATime,surv.event = data$STIA %>% as.numeric,method="noether",na.rm = T)
c_index_1 <- c(str_c(c(c_index$c.index %>% round(3),"(",c_index$lower %>% round(3),"-",c_index$upper %>% round(3),")"),collapse = ""))        
}) 
cindex_with_ci

cindex <- fit_ti_rocs %>% lapply(function(x){c_index <- concordance.index(predict(x,data),
                                                                          surv.time = data$STIATime,surv.event = data$STIA %>% as.numeric,method="noether",na.rm = T)
#cindex_with_ci <- cindex_with_ci %>% data.frame()

#c_index_1 <- c_index$c.index %>% round(3)    
})

## C-statistics 比較,比較的是列表，不是單個值
cind_comp_p <- c()
i = 1
while (i < length(cindex)) {
  p <- cindex.comp(cindex[[i+1]],cindex[[1]])
  p <- p[[1]] 
  cind_comp_p <- append(cind_comp_p,p)
  i = i+1
}
cind_comp_p

七、計算并比較cox模型IDI和NRI

##-----------------------------計算cox模型IDI和NRI--------------------------------------------
library(survIDINRI)
#定義時間節(jié)點
#t0=365*5
t0=60
## 所有列轉成數(shù)字
data_ti_roc <- apply(data_ti_roc,2,function(x)as.numeric(x)) %>% data.frame()

#data_ti_roc <- data_ti_roc[1:100,]       ## 先用小樣本實驗，計算時間很長

## 提取結局列
outcomes = as.matrix(subset(data_ti_roc,select=c(STIATime,STIA)))

##提取變量列，轉成矩陣

#data_ti_roc$time = as.integer(data_ti_roc$STIATime)
#data_ti_roc$status = as.numeric(data_ti_roc$STIA)
head(data_ti_roc)
model_1 <- as.matrix(subset(data_ti_roc,select=c("CHA2DS2VASC"))) 
model_2 <- as.matrix(subset(data_ti_roc,select=c("CHA2DS2VASC","PTFV15000"))) 
model_3 <- as.matrix(subset(data_ti_roc,select=c(CHA2DS2VASC,LAE))) 
model_4 <- as.matrix(subset(data_ti_roc,select=c(CHA2DS2VASC,LAE,PTFV15000))) 
head(model_1)
head(model_2)
head(model_3)
head(model_4)
head(data_ti_roc)

##開始對比

x2<-IDI.INF(outcomes,model_1, model_2, t0, npert=1000)
IDI.INF.OUT(x2)
x3<-IDI.INF(outcomes,model_1, model_3, t0, npert=1000)
IDI.INF.OUT(x3)
x4<-IDI.INF(outcomes,model_1, model_4, t0, npert=1000)
IDI.INF.OUT(x4)
#只能手抄結果
#M1表示IDI
#M2表示NRI
#M3表示中位數(shù)差異

八、繪制COX校正曲線

##----------------------------制作cox校正曲線-----------------------------
library(rms)

dd<-datadist(data_ti_roc) #設置工作環(huán)境變量，將數(shù)據(jù)整合
options(datadist='dd') #設置工作環(huán)境變量，將數(shù)據(jù)整合
##
##time.in 和 u 要是一樣的，都是要評價的時間節(jié)點
time_inc <- 60
coxm_1 <- cph(formula_ti_roc_4,data=data_ti_roc,surv=T,x=T,y=T,time.inc = time_inc)
###繪制cox回歸生存概率的nomogram圖
## 構建Nomo圖的對象只能是rms保重d額cph()函數(shù)
surv <- Survival(coxm_1) # 構建生存概率函數(shù)
nom <- nomogram(coxm_1,fun=list(#function(x)surv(70,1-x),
                                function(x)surv(60,1-x),
                                function(x)surv(36,1-x)),  ##算出不同時間節(jié)點生存率值
                ##算出不同時間節(jié)點生存率值，顯示在列線圖上
                funlabel = c("60","36"),
                lp=F
                #fun.at=c('0.9','0.85','0.80','0.70','0.6','0.5','0.4','0.3','0.2','0.1')
)
plot(nom)

##time.in 和 u 要是一樣的，都是要評價的時間節(jié)點
time_inc <- 60
coxm_1 <- cph(formula_ti_roc_4,data=data_ti_roc,surv=T,x=T,y=T,time.inc = time_inc)
m = (nrow(data_ti_roc)/5) %>% ceiling()   ## m 約等于樣本數(shù)5分之一
m
cal_1<-calibrate(coxm_1,u=time_inc,cmethod='KM',m= m,B=1000)

pdf("fig5a.pdf",width = 3.9,height = 4.5)

plot(cal_1,lwd=2,lty=1, ##設置線條形狀和尺寸
     errbar.col=c(rgb(0,118,192,maxColorValue = 255)), ##設置一個顏色
     xlab='Nomogram-Predicted Probability of 1-year DFS',#便簽
     ylab='Actual 1-year DFS(proportion)',#標簽
     col=c(rgb(192,98,83,maxColorValue = 255)),#設置一個顏色
     xlim = c(0,1),ylim = c(0,1)) ##x軸和y軸范圍
dev.off()

table(data$STIA)

九、制作cox 決策曲線

##-----------------------------制作cox 決策曲線-----------------------------------------
library(ggDCA)
if (!require(ggDCA)) {
  devtools::install_github("yikeshu0611/ggDCA")
}

dca_1 <- dca(fit_ti_roc_3)
dca_4 <- dca(fit_ti_roc_4)
dca <- dca(fit_ti_roc_1,fit_ti_roc_4,times = 60,model.names = c("model1","model4"))
dca_1 <- ggplot(dca,       
                model.names="模型1",
                linetype =F, #線型
                lwd = 0.6)+
  #ylim(c(-0.02,0.05))+
  theme(legend.position = c(0.7,0.8),
        text = element_text(size = 8)
        )+
  scale_color_lancet()
dca_1

ggsave("time_dca.pdf",dca_1,width = 8,height = 8.5,units = "cm")