library(pROC)  #對(duì)于AUC值，采用pROC包計(jì)算
library(rms)


####################
#訓(xùn)練集Vdata2
#LVEF劃分＜50；；≥50   √
Vdata2$LVEF_F<- cut(Vdata2$LVEF, breaks = c(-Inf, 49.9, Inf),
                    labels = c("0", "1"), include.lowest = TRUE)

summary(Vdata2)
####################



#lrm得C指數(shù)
fit3 <- lrm(MACE ~ Age +Killip +Hypertention+Hhcy +HF +Sent+VD3LM
            +LDLC+LVEF_F,x = TRUE,y = TRUE,data = Vdata2)

#glm得ROC曲線

#
fit4 <- glm(MACE ~ Age +Killip +Hypertention+Hhcy +HF +Sent+VD3LM
            +LDLC+LVEF_F,data = Vdata2,family = "binomial")

logistic.display(fit4) 

pre <- predict(fit4,type='response')
plot.roc(Vdata2$MACE, pre,
         main="ROC Curve", percent=TRUE,
         print.auc=TRUE,
         ci=TRUE, of="thresholds",
         thresholds="best",
         print.thres="best",
         legacy.axes=TRUE)


#計(jì)算AUC的置信區(qū)間
rocplot1 <- roc(Vdata2$MACE,pre)
ci.auc(rocplot1)


##############備選代碼作圖

######ROC 基于plotroc包
model.train<-glm(formula = MACE ~ Age +Killip +Hypertention+Hhcy +HF +Sent+VD3LM
                 +LDLC+LVEF_F, family = binomial(), 
                 data = Vdata2)


Vdata2$p_prediction <- predict(model.train, type="response")
roc.train <- roc(Vdata2$MACE, Vdata2$p_prediction)
roc.train
plot(roc.train)

#或用ggplot包
library(ggplot2)

rocplot <- ggplot(Vdata2, aes(d = MACE  , m = p_prediction)) + geom_roc(n.cuts = 0) + 
  style_roc(theme = theme_bw, xlab = "1-Specificity", ylab = "Sensitivity") + 
  geom_abline(intercept = 0, slope = 1) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
rocplot


######################



#校準(zhǔn)曲線
fit3 <- lrm(Vdata2$MACE ~ pre,x = TRUE,y = TRUE)
cal2 <- calibrate(fit3,  method = "boot", B = 1000)
plot(cal2, xlab = "Predicted Probability", ylab = "Actual Probability",main = "Calibration Curve")



# calibration plot

val.prob(p = fitted(model.train), y = Vdata2$MACE,logistic.cal=F)

###################################################################33





#驗(yàn)證集Vdata1

#LVEF劃分＜50；；≥50   √
Vdata1$LVEF_F<- cut(Vdata1$LVEF, breaks = c(-Inf, 49.9, Inf),
                    labels = c("0", "1"), include.lowest = TRUE)

summary(Vdata1)
##################################################3

library(pROC)  #對(duì)于AUC值，采用pROC包計(jì)算

#lrm得C指數(shù)
fit5 <- lrm(MACE ~ Age +Killip +Hypertention+Hhcy +HF +Sent+VD3LM
            +LDLC+LVEF_F,x = TRUE,y = TRUE,data = Vdata1)



#glm得ROC曲線

#
fit6 <-  glm(MACE ~ Age +Killip +Hypertention+Hhcy +HF +Sent+VD3LM
             +LDLC+LVEF_F,data = Vdata1,family = "binomial")



pre <- predict(fit6,type='response')
plot.roc(Vdata1$MACE, pre,
         main="ROC Curve", percent=TRUE,
         print.auc=TRUE,
         ci=TRUE, of="thresholds",
         thresholds="best",
         print.thres="best",
         legacy.axes=TRUE)


#計(jì)算AUC的置信區(qū)間
rocplot1 <- roc(Vdata1$MACE,pre)
ci.auc(rocplot1)




#備選作圖######ROC 基于plotroc包
model.train<-glm(formula = MACE ~ Age +Killip +Hypertention+Hhcy +HF +Sent+VD3LM
                 +LDLC+LVEF_F+InhP2Y12, family = binomial(), 
                 data = Vdata1)


Vdata1$p_prediction <- predict(model.train, type="response")
roc.train <- roc(Vdata1$MACE, Vdata1$p_prediction)
roc.train
plot(roc.train)

#或用ggplot包
library(ggplot2)

rocplot <- ggplot(Vdata1, aes(d = MACE  , m = p_prediction)) + geom_roc(n.cuts = 0) + 
  style_roc(theme = theme_bw, xlab = "1-Specificity", ylab = "Sensitivity") + 
  geom_abline(intercept = 0, slope = 1) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
rocplot


######################



#校準(zhǔn)曲線
fit5 <- lrm(MACE ~ Age +Killip +Hypertention+Hhcy +HF +Sent+VD3LM
            +LDLC+LVEF_F+InhP2Y12,x = TRUE,y = TRUE,data = Vdata1)

fit5 <- lrm(Vdata1$MACE ~ pre,x = TRUE,y = TRUE)
cal2 <- calibrate(fit5,  method = "boot", B = 1000)
plot(cal2, xlab = "Predicted Probability", ylab = "Actual Probability",main = "Calibration Curve")