#install.packages("foreign")
#install.packages("rms")
#install.packages("pROC")
#install.packages("rmda")
#install.packages("nricens")
library(foreign)
library(rms)   #构建Logstic模型
library(pROC)   #绘制ROC曲线
library(rmda)   #绘制DCA曲线
library(nricens)   #计算NRI值

## setwd("C:\Users\000\Desktop\14_Logstic")   #设置工作目录
# rt <- read.table("C:\\Users\\yulei\\Desktop\\0120.txt",header=T,sep="\t")   #读取数据
rt <- read.csv("C:\\Users\\yulei\\Desktop\\0306.csv",header=TRUE)
## LIHC <- read.table(file = 'TCGA-LIHC-nomogram.tsv', sep = '\t', header = TRUE)
head(rt)


## 关键步骤：按照nomogram要求“打包”数据，可以输入??datadist查看详细说明
dd=datadist(rt)
options(datadist="dd") 

fit <- glm(group ~ Age + RBC + Hb + Hct + Fibrous.calcification, data=rt, 
           family = binomial()) 

#利用lrm(）函数对模型进行拟合 

summary(fit)   #查看模型拟合结果

Cindex <- rcorrcens(group~predict(fit), data = rt)
Cindex

cal <- calibrate(gfit, method="boot", B=1000)
plot(cal,
     xlab="Nomogram-predicted probability of nonadherence",
     ylab="Actual diagnosed nonadherence (proportion)",
     sub=F)


# 全局结果roc曲线绘制
gfit <- roc(group~predict(fit), data = rt)
plot(gfit,
     print.auc=TRUE,   #输出AUC值
     print.thres=FALSE,   #输出cut-off值
     main = "ROC CURVE",   #设置图形的标题
     col= "red",   #曲线颜色
     print.thres.col="black",   #cut-off值字体的颜色
     identity.col="lightblue",   #对角线颜色
     identity.lty=1,identity.lwd=1,
     legacy.axes=TRUE)


# 全局加单因素roc曲线绘制
gfit <- roc(group~predict(fit), data = rt,cause=1)
roc_age <-roc(group~Age,data = rt)
roc_rbc <-roc(group~RBC,data = rt)
roc_hb <-roc(group~Hb,data = rt)
roc_hct <-roc(group~Hct,data = rt)
roc_fc <-roc(group~Fibrous.calcification,data = rt)

par(mar= c(5,5,1,1),cex.lab=1,cex.axis= 1) #设置图形边界

# plot.roc(gfit,col="red",precision=seq(0.1,0.9, by=0.1) )
# axis(1,seq(0, 1, 0.2), seq(0, 2, 0.2))

plot(gfit,
     print.auc=TRUE,   #输出AUC值
     print.thres=FALSE,   #输出cut-off值
     main = "ROC CURVE",   #设置图形的标题
     col= "red",   #曲线颜色
     print.thres.col="grey",   #cut-off值字体的颜色
     identity.col="pink",   #对角线颜色
     identity.lty=1,identity.lwd=1,
     legacy.axes=TRUE)
lines.roc(roc_age, col= "purple")
#lines.roc(roc_rbc, col= "blue")
#lines.roc(roc_hb, col= "black")
#lines.roc(roc_hct, col= "green")
lines.roc(roc_fc, col= "green")

legend(x = 0.35, y = 0.25,
       legend = c("nomogram","age", "Fibrous_calcification"),
       fill = c("red", "purple", "green"),
       cex = 0.8)
#plot(gfit, identify

summary(gfit)