########## 1. Install and load required libraries ##########
#install.packages(c("caret", "readxl", "pROC", "RcmdrPlugin.EZR", "ggplot2")) ##This needs to bedone only once.
library(caret, readxl, pROC)
library(RcmdrPlugin.EZR)
require("ggplot2")
setwd("~/Downloads")

########## 2. Load Excel Spreadsheet and recode the target variable ##########

Dataset <- read_excel("ImagingFeatures.xlsx", sheet = "Export")
Dataset$Tumor <- make.names(Dataset$Tumor)

########## 3. Determine the optimal hyperparameter combination ##########

TunedGrid <- expand.grid(size=(1:7),
                         decay=c(0.0001,0.0005,0.001,0.005,0.01),
                         bag=F)
set.seed(1234)
Testmodel= train(as.factor(Tumor) ~ AUC+PE+SUVmax, 
                 data = Dataset, 
                 method = "avNNet", 
                 preProc=c("BoxCox"),
                 tuneGrid = TunedGrid, 
                 trControl = trainControl(method="boot",classProbs = TRUE, savePredictions=F))
Testmodel$bestTune

########## 4. Calculate the final model ##########

TunedGrid <- expand.grid(size=Testmodel$bestTune$size,
                         decay=Testmodel$bestTune$decay,
                         bag=F)
set.seed(1234)
NeuralNetModel= train(as.factor(Tumor) ~ AUC+PE+SUVmax, 
                      data = Dataset, 
                      method = "avNNet", 
                      preProc=c("BoxCox"),
                      tuneGrid = TunedGrid, 
                      trControl = trainControl(method="LOOCV", classProbs = TRUE, savePredictions=T))
NeuralNetModel

########## 5. Calculate standard parameters of diagnostic accuracy ##########

epitable <- confusionMatrix(NeuralNetModel$pred[,1], NeuralNetModel$pred[,2], positive="X1") 
epi.tests(rev(epitable$table), conf.level = 0.95)

########## 6. Plot the ROC curve ##########

Predictions <- as.data.frame(NeuralNetModel$pred)
roc.obj <- pROC::roc(obs ~ X1, data=Predictions)
plot(roc.obj, print.auc=T)

########## 7. Compare the final model's ROC curve with the ROC curves of its constituent parameters  ##########
Dataset$NNPrediction <- Predictions$X1

ROC.NN <- pROC::roc(Tumor ~ NNPrediction, data=Dataset, na.rm=TRUE, percent=FALSE, direction='auto',                    partial.auc=FALSE, partial.auc.focus='specificity', partial.auc.correct=FALSE, auc=TRUE, plot=FALSE,                     ci=TRUE, of='auc', conf.level=0.95, ci.method='bootstrap', boot.n=2000, boot.stratified=FALSE)
ROC.AUC <- pROC::roc(Tumor ~ AUC, data=Dataset, na.rm=TRUE, percent=FALSE, direction='auto', partial.auc=FALSE, partial.auc.focus='specificity', partial.auc.correct=FALSE, auc=TRUE, plot=FALSE, ci=TRUE, of='auc', conf.level=0.95, ci.method='bootstrap', boot.n=2000, boot.stratified=FALSE)
ROC.PE <- pROC::roc(Tumor ~ PE, data=Dataset, na.rm=TRUE, percent=FALSE, direction='auto', partial.auc=FALSE, partial.auc.focus='specificity', partial.auc.correct=FALSE, auc=TRUE, plot=FALSE, ci=TRUE, of='auc', conf.level=0.95, ci.method='bootstrap', boot.n=2000, boot.stratified=FALSE)
ROC.SUVmax <- pROC::roc(Tumor ~ SUVmax, data=Dataset, na.rm=TRUE, percent=FALSE, direction='auto', partial.auc=FALSE, partial.auc.focus='specificity', partial.auc.correct=FALSE, auc=TRUE, plot=FALSE, ci=TRUE, of='auc', conf.level=0.95, ci.method='bootstrap', boot.n=2000, boot.stratified=FALSE)

roc.test(ROC.NN,ROC.AUC,paired=TRUE,alternative='two.sided', boot.n=2000)
roc.test(ROC.NN,ROC.PE, paired=TRUE,alternative='two.sided', boot.n=2000)
roc.test(ROC.NN,ROC.SUVmax,paired=TRUE,alternative='two.sided', boot.n=2000)