Script_MLvsLR.Rmd

---
title: "PONES-21-39032_Script"
author: "Sara Dominguez"
date: '2022-02-28'
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```



```{r remove}
rm(list=ls())
graphics.off()
```


    # L I B A R I E S
```{r}
library(missForest)
library(caret)
library(caretEnsemble)
library(dplyr)
library(plyr)
library(MASS)
library(plotROC)
library(MLeval)
library(extrafont)
loadfonts()
library(ggpubr)
library(readr)
library(reshape2)
library(pROC)
library(multcomp)
```



    # D A T A
    
Data is not publicly available because is protected by European GDPR (data pseudoanonymized). However, can be formally shared under a formal application and research proposal after PENTA Foundation acceptance. Please send your proposal to sara.dominguez.r@gmail.com
    
    MAY 2020
```{r}
#load("path/mydata_dp12.rda")
```
    
    Code the data
```{r}
mydata_dp12<-as.data.frame(mydata_dp12)
mydata_dp12$birth_sex<-recode_factor(mydata_dp12$birth_sex, "1"="Female", "2"="Male")
mydata_dp12$dp12<-recode_factor(mydata_dp12$dp12, "1"="Yes", "0"="No")
mydata_dp12$WAZ<-as.numeric(mydata_dp12$WAZ$waz)
mydata_dp12$preterm_birth<-as.factor(mydata_dp12$preterm_birth)
mydata_dp12$children_baseline_ART<-as.factor(mydata_dp12$children_baseline_ART)
mydata_dp12$mother_SLEHI<-as.factor(mydata_dp12$mother_SLEHI)

mydata_dp12<-mydata_dp12[,c(1:12,15)]
```
    

    # D A T A  S P L I T
```{r}

#train and test
#training
train.rows<- createDataPartition(y= mydata_dp12$dp12, p=0.7, list = FALSE)
train.data12<- mydata_dp12[train.rows,] # 70% data goes in here
summary(factor(train.data12$dp12))

#testing
test.data12<- mydata_dp12[-train.rows,] # 30% data goes in here
summary(factor(test.data12$dp12))

setwd("C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Data/")
save(train.data12,file = paste(paste("train.data_dp12",Sys.Date(),sep = "_"),".Rda",sep=""))
save(test.data12,file = paste(paste("test.data_dp12",Sys.Date(),sep = "_"),".Rda",sep=""))

```
    


    # I M P U T E
```{r}
#DP 12
colnames(train.data12)
train.data12_imp <- missForest(train.data12[,c(2:12)])
train.data12_imp$OOBerror
train.data12_i<- cbind(PATIENT=train.data12$record_id,train.data12_imp$ximp, dp12=train.data12$dp12)


#na omit
train.data12_completecase<-na.omit(train.data12)
summary(factor(train.data12_completecase$dp12))

save(train.data12_i, file = paste(paste("train.data12_imputed",Sys.Date(),sep = "_"),".Rda",sep=""))
save(train.data12_completecase, file = paste(paste("train.data12_cc",Sys.Date(),sep = "_"),".Rda",sep=""))


test.data12_imp <- missForest(test.data12[,c(2:12)])
test.data12_imp$OOBerror
test.data12_i<- cbind(PATIENT=test.data12$record_id,test.data12_imp$ximp, dp12=test.data12$dp12)

#na omit
test.data12_completecase<-na.omit(test.data12)
summary(factor(test.data12_completecase$dp12))

save(test.data12_i, file = paste(paste("test.data12_imputed",Sys.Date(),sep = "_"),".Rda",sep=""))
save(test.data12_completecase, file = paste(paste("test.data12_cc",Sys.Date(),sep = "_"),".Rda",sep=""))
```
    
    
    
    # T R A I N  C O N T R O L

# Train control crossvalidation

```{r}
#crossval
set.seed(12345)

ctrl <- trainControl(method = "repeatedcv", repeats = 5,
                     classProbs = TRUE,
                     summaryFunction = twoClassSummary,
                     ## new option here:
                     sampling = "down", savePredictions = T)
```



    
#-------------------------------------------------------------------------------------------
    # DP12: 12-months Death and Progression
#-------------------------------------------------------------------------------------------
 

    # Model 1: Logistic Regression

```{r}
mod<-glm(dp12 ~ .,data = train.data12_i[,c(2:13)], family=binomial)
stepAIC(mod)


glmFit <- train(dp12 ~ .,data = train.data12_i[,c(2:13)],method = "glm", family="binomial",trControl=ctrl)

glmClasses <- predict(glmFit, newdata = test.data12_i)
glmProbs <- predict(glmFit, newdata = test.data12_i, type = "prob")

confusionMatrix(data = glmClasses , test.data12_i$dp12)
saveRDS(glmFit, "C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Models/glmFit.rds")


#system time
system.time(glmFit <- train(dp12 ~ .,data = train.data12_i[,c(2:13)],method = "glm", family="binomial",trControl=ctrl)) #1.48

time_glm<-c(1.36, 1.35,1.35,1.31,1.39,1.31,1.39,1.35,1.33,1.32)
summary(time_glm)



#AUC
result.roc <- roc(test.data12_i$dp12, glmProbs$Yes)

#Variable importance
gbmImp <- varImp(glmFit, scale = FALSE)
plot(gbmImp)
```


        # Model 2: Random Forest
```{r}
rfFit2 <- train(dp12 ~ .,data = train.data12_i[,c(2:13)],method = "rf",trControl=ctrl, tuneLength=10)
rdaGrid = data.frame(mtry = 12)
rfFit <- train(dp12 ~ .,data = train.data12_i[,c(2:13)],method = "rf",trControl=ctrl, tuneGrid = rdaGrid)

rfClasses <- predict(rfFit, newdata = test.data12_i)
rfProbs <- predict(rfFit, newdata = test.data12_i, type = "prob")

confusionMatrix(data = rfClasses , test.data12_i$dp12) #sensitiviy 0.7 #specificity 0.8 
saveRDS(rfFit, "C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Models/rfFit.rds")

vari<-varImp(rfFit$finalModel)
varImpPlot(rfFit$finalModel)


system.time(rfFit <- train(dp12 ~ .,data = train.data12_i[,c(2:13)],method = "rf",trControl=ctrl, tuneGrid = rdaGrid)) #1.89

time_rf<-c(1.95, 1.98, 1.92, 2, 2, 2.03, 2.0, 1.98, 1.96, 2.0)
summary(time_rf)

#Variable importance
rfImp <- varImp(rfFit, scale = FALSE)
plot(rfImp)
```
        

        # Model 3: Support Vector Machine
```{r}
svmFit <- train(dp12 ~ .,data = train.data12_i[,c(2:13)],method = "svmRadial",trControl=ctrl, tuneLength=10)
rdaGrid = data.frame(C = 8, sigma= 4.694759e-11)
svmFit <- train(dp12 ~ .,data = train.data12_i[,c(2:13)],method = "svmRadial",trControl=ctrl, tuneGrid = rdaGrid)

svmClasses <- predict(svmFit, newdata = test.data12_i)
svmProbs <- predict(svmFit, newdata = test.data12_i, type = "prob")

confusionMatrix(data = svmClasses , test.data12_i$dp12) 

saveRDS(svmFit, "C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Models/svmFit.rds")


system.time(svmFit <- train(dp12 ~ .,data = train.data12_i[,c(2:13)],method = "svmRadial",trControl=ctrl, tuneGrid = rdaGrid) #1.63
)

time_svm<-c(2.17,1.65, 1.68, 1.7, 1.68, 1.71, 1.7, 1.68,1.7,1.7)
summary(time_svm)

#Variable importance
svmImp <- varImp(svmFit, scale = FALSE)
plot(svmImp)
```
        

          # Model 4: Naive Bayes
```{r}
nbFit <- train(dp12 ~ .,data = train.data12_i[,c(2:7,9:13)],method = "nb")
rdaGrid = data.frame(fL=0, usekernel=TRUE, adjust=1)
nbFit <- train(dp12 ~ .,data = train.data12_i[,c(2:7,9:13)],method = "nb",trControl=ctrl, tuneGrid = rdaGrid)

nbClasses <- predict(nbFit, newdata = test.data12_i)
nbProbs <- predict(nbFit, newdata = test.data12_i, type = "prob")

confusionMatrix(data = nbClasses , test.data12_i$dp12) #sensitiviy 0.8 #specificity 0.7


saveRDS(nbFit, "C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Models/nbFit.rds")


system.time(nbFit <- train(dp12 ~ .,data = train.data12_i[,c(2:7,9:13)],method = "nb",trControl=ctrl, tuneGrid = rdaGrid)) #2.47

time_nb<-c(2.5,2.69,2.58,2.61, 2.83, 2.52, 2.61, 2.56,2.64, 2.61)
summary(time_nb)

#Variable importance
nbImp <- varImp(nbFit, scale = FALSE)
plot(nbImp)
```
        
        # Model 5: KNN
```{r}
knnFit <- train(dp12 ~ .,data = train.data12_i[,c(2:13)],method = "knn",trControl=ctrl, tuneLength=10)
rdaGrid = data.frame(k=5)
knnFit <- train(dp12 ~ .,data = train.data12_i[,c(2:13)],method = "knn",trControl=ctrl, tuneGrid = rdaGrid)

knnClasses <- predict(knnFit, newdata = test.data12_i)
knnProbs <- predict(knnFit, newdata = test.data12_i, type = "prob")

confusionMatrix(data = knnClasses , test.data12_i$dp12) 

saveRDS(knnFit, "C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Models/knnFit.rds")


system.time(knnFit <- train(dp12 ~ .,data = train.data12_i[,c(2:13)],method = "knn",trControl=ctrl, tuneGrid = rdaGrid)) #0.95

time_knn<-c(0.99,1.05,1.01,1.05,1.06,1.05, 1.07,1.04,1.05,1)
summary(time_knn)

#Variable importance
knnImp <- varImp(knnFit, scale = FALSE)
plot(knnImp)
```
        
    # Model 6: ANN
```{r}
annFit <- train(dp12 ~ .,data = train.data12_i[,c(2:13)],method = "nnet",trControl=ctrl, tuneLength=10)
rdaGrid = data.frame(size=11, decay=0.1)
annFit <- train(dp12 ~ .,data = train.data12_i[,c(2:13)],method = "nnet",trControl=ctrl, tuneGrid = rdaGrid)

annClasses <- predict(annFit, newdata = test.data12_i)
annProbs <- predict(annFit, newdata = test.data12_i, type = "prob")

confusionMatrix(data = annClasses , test.data12_i$dp12) 

saveRDS(annFit, "C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Models/annFit.rds")

system.time(annFit <- train(dp12 ~ .,data = train.data12_i[,c(2:13)],method = "nnet",trControl=ctrl, tuneGrid = rdaGrid) #2.3
)

time_ann<-c(2.58, 2.45, 2.5,2.56,2.54, 2.5,2.5, 2.5, 2.5, 2.43)
summary(time_ann)

#Variable importance
annImp <- varImp(annFit, scale = FALSE)
plot(annImp)
```
   
   
       # Model 7: GLMNET
```{r}
glmnetFit <- train(dp12 ~ .,data = train.data12_i[,c(2:13)],method = "glmnet",trControl=ctrl, tuneLength=10)
rdaGrid = data.frame(alpha=0.8, lambda=0.210781)
glmnetFit <- train(dp12 ~ .,data = train.data12_i[,c(2:13)],method = "glmnet",trControl=ctrl, tuneGrid = rdaGrid)

glmnetClasses <- predict(glmnetFit , newdata = test.data12_i)
glmnetProbs <- predict(glmnetFit , newdata = test.data12_i, type = "prob")

confusionMatrix(data = glmnetClasses , test.data12_i$dp12) #sensitiviy 0.8 #specificity 0.7

saveRDS(glmnetFit, "C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Models/glmnetFit.rds")

system.time(glmnetFit <- train(dp12 ~ .,data = train.data12_i[,c(2:13)],method = "glmnet",trControl=ctrl, tuneGrid = rdaGrid)) #1.81

time_glmnet<-c(1.95,1.75,1.79, 1.83, 1.87, 1.78, 1.81, 1.8,1.8,1.86)
summary(time_glmnet)

#Variable importance
glmnetImp <- varImp(glmnetFit, scale = FALSE)
plot(glmnetImp)

```
   
  
   
   
   
       # M o d e l s  R O C
```{r}
p1<-roc(glmFit$pred$obs,glmFit$pred$Yes,auc = T,
            smooth = F,
            # arguments for ci
            ci=TRUE, ci.alpha=0.9, stratified=FALSE,
            # arguments for plot
            plot=TRUE, auc.polygon=FALSE, grid=TRUE,
            print.auc=TRUE, show.thres=F, percent=TRUE, col="brown4",
            asp = 0, print.thres=F, main = "Logistic regression", legacy.axes=T)
glmroc<-ggroc(p1, colour="aquamarine3", size=1.5) + theme(panel.background = element_blank(), panel.border = element_rect(fill=NA, colour="black"), text=element_text(family="Segoe UI Historic", size=12)) + annotate(geom="text", x=30, y=25, label="AUC: 65.3% (58.9%-71.8%)", family="Segoe UI Historic") + ggtitle("Logistic")

#optimisim glm
glmProbs_train <- predict(glmFit, newdata = train.data12_i, type = "prob")
AUC_apparent<-roc(train.data12_i$dp12, glmProbs_train$Yes, ci=T)
AUC_apparent<-AUC_apparent$auc
AUC_validated<-roc(test.data12_i$dp12, glmProbs$Yes, ci=T) #0.68 (0.4648-0.8908)
AUC_validated<-AUC_validated$auc
optimisim_glm<-AUC_apparent - AUC_validated


p2<-roc(rfFit$pred$obs,rfFit$pred$Yes,auc = T,
            smooth = F,
            # arguments for ci
            ci=TRUE, ci.alpha=0.9, stratified=FALSE,
            # arguments for plot
            plot=TRUE, auc.polygon=FALSE, grid=TRUE,
            print.auc=TRUE, show.thres=F, percent=TRUE, col="brown4",
            asp = 0, print.thres=F, main = "Random Forest", legacy.axes=T)
rfroc<-ggroc(p2, colour="salmon", size=1.5) + theme(panel.background = element_blank(), panel.border = element_rect(fill=NA, colour="black"), text=element_text(family="Segoe UI Historic", size=12)) + annotate(geom="text", x=40, y=50, label="AUC: 73.2% (67.2%-79.1%)", family="Segoe UI Historic") + ggtitle("Random Forest")

#optimisim rf
rfProbs_train <- predict(rfFit, newdata = train.data12_i, type = "prob")
AUC_apparent<-roc(train.data12_i$dp12, rfProbs_train$Yes, ci=T)
AUC_apparent<-AUC_apparent$auc
AUC_validated<-roc(test.data12_i$dp12, rfProbs$Yes, ci=T) #0.68 (0.4648-0.8908)
AUC_validated<-AUC_validated$auc
optimisim_rf<-AUC_apparent - AUC_validated



p3<-roc(svmFit$pred$obs,svmFit$pred$Yes,auc = T,
            smooth = F,
            # arguments for ci
            ci=TRUE, ci.alpha=0.9, stratified=FALSE,
            # arguments for plot
            plot=TRUE, auc.polygon=FALSE, grid=TRUE,
            print.auc=TRUE, show.thres=F, percent=TRUE, col="brown4",
            asp = 0, print.thres=F, main = "SVM", legacy.axes=T)
svmroc<-ggroc(p3, colour="purple", size=1.5) + theme(panel.background = element_blank(), panel.border = element_rect(fill=NA, colour="black"), text=element_text(family="Segoe UI Historic", size=12)) + annotate(geom="text", x=20, y=25, label="AUC: 46.0% (39.9%-52.2%)", family="Segoe UI Historic")  + ggtitle("Support Vector Machine")

#optimisim svm
svmProbs_train <- predict(svmFit, newdata = train.data12_i, type = "prob")
AUC_apparent<-roc(train.data12_i$dp12, svmProbs_train$Yes, ci=T)
AUC_apparent<-AUC_apparent$auc
AUC_validated<-roc(test.data12_i$dp12, svmProbs$Yes, ci=T) #0.68 (0.4648-0.8908)
AUC_validated<-AUC_validated$auc
optimisim_svm<-AUC_apparent - AUC_validated


p4<-roc(nbFit$pred$obs,nbFit$pred$Yes,auc = T,
            smooth = F,
            # arguments for ci
            ci=TRUE, ci.alpha=0.9, stratified=FALSE,
            # arguments for plot
            plot=TRUE, auc.polygon=FALSE, grid=TRUE,
            print.auc=TRUE, show.thres=F, percent=TRUE, col="brown4",
            asp = 0, print.thres=F, main = "SVM", legacy.axes=T)
nbroc<-ggroc(p4, colour="cornflowerblue", size=1.5) + theme(panel.background = element_blank(), panel.border = element_rect(fill=NA, colour="black"), text=element_text(family="Segoe UI Historic", size=12)) + annotate(geom="text", x=30, y=50, label="AUC: 65.1% (58.7%-71.6%)", family="Segoe UI Historic")  + ggtitle("Naïve Bayes")

#optimisim nb
nbProbs_train <- predict(nbFit, newdata = train.data12_i, type = "prob")
AUC_apparent<-roc(train.data12_i$dp12, nbProbs_train$Yes, ci=T)
AUC_apparent<-AUC_apparent$auc
AUC_validated<-roc(test.data12_i$dp12, nbProbs$Yes, ci=T) #0.68 (0.4648-0.8908)
AUC_validated<-AUC_validated$auc
optimisim_nb<-AUC_apparent - AUC_validated



p5<-roc(knnFit$pred$obs,knnFit$pred$Yes,auc = T,
            smooth = F,
            # arguments for ci
            ci=TRUE, ci.alpha=0.9, stratified=FALSE,
            # arguments for plot
            plot=TRUE, auc.polygon=FALSE, grid=TRUE,
            print.auc=TRUE, show.thres=F, percent=TRUE, col="brown4",
            asp = 0, print.thres=F, main = "SVM", legacy.axes=T)
knnroc<-ggroc(p5, colour="darkgoldenrod2", size=1.5) + theme(panel.background = element_blank(), panel.border = element_rect(fill=NA, colour="black"), text=element_text(family="Segoe UI Historic", size=12)) + annotate(geom="text", x=40, y=50, label="AUC: 69.6% (63.2%-76.1%)", family="Segoe UI Historic")  + ggtitle("K-nearest neighbour")

#optimisim knn
knnProbs_train <- predict(knnFit, newdata = train.data12_i, type = "prob")
AUC_apparent<-roc(train.data12_i$dp12, knnProbs_train$Yes, ci=T)
AUC_apparent<-AUC_apparent$auc
AUC_validated<-roc(test.data12_i$dp12, knnProbs$Yes, ci=T) #0.68 (0.4648-0.8908)
AUC_validated<-AUC_validated$auc
optimisim_knn<-AUC_apparent - AUC_validated



p6<-roc(annFit$pred$obs,annFit$pred$Yes,auc = T,
            smooth = F,
            # arguments for ci
            ci=TRUE, ci.alpha=0.9, stratified=FALSE,
            # arguments for plot
            plot=TRUE, auc.polygon=FALSE, grid=TRUE,
            print.auc=TRUE, show.thres=F, percent=TRUE, col="brown4",
            asp = 0, print.thres=F, main = "SVM", legacy.axes=T)
annroc<-ggroc(p6, colour="hotpink", size=1.5) + theme(panel.background = element_blank(), panel.border = element_rect(fill=NA, colour="black"), text=element_text(family="Segoe UI Historic", size=12)) + annotate(geom="text", x=30, y=30, label="AUC: 52.9% (46.4%-59.4%)", family="Segoe UI Historic")  + ggtitle("Neural Network")

#optimisim ann
annProbs_train <- predict(annFit, newdata = train.data12_i, type = "prob")
AUC_apparent<-roc(train.data12_i$dp12, annProbs_train$Yes, ci=T)
AUC_apparent<-AUC_apparent$auc
AUC_validated<-roc(test.data12_i$dp12, annProbs$Yes, ci=T) #0.68 (0.4648-0.8908)
AUC_validated<-AUC_validated$auc
optimisim_ann<-AUC_apparent - AUC_validated

ggarrange(glmroc, rfroc, svmroc, nbroc, knnroc, annroc, nrow=3, ncol=2)

roc.test(p1, p2)

```
    
    
    # P e r f o r m a n c e  P l o t
```{r}
Models_Performance <- read_delim("C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Models_Performance_24-05-2020.csv", 
    ";", escape_double = FALSE, locale = locale(decimal_mark = ","), 
    trim_ws = TRUE)

Models_Performance<-Models_Performance[Models_Performance$Metric!="Kappa", ]

Models_Performance$Model<-factor(Models_Performance$Model, levels=c("Logistic","Glmnet", "Random Forest", "Support Vector Machine", "Naive Bayes", "K-Nearest Neighbor", "Artificial Neural Networks"))

Models_Performance$Metric<-factor(Models_Performance$Metric, levels=rev(c("Accuracy", "Sensitivity", "Specificity", "PPV", "NPV")))

p = ggplot(data=Models_Performance,
    aes(x = Metric,y = Value, ymin = lower, ymax = upper ))+
    geom_pointrange(aes(col=Model), position = position_dodge(0.6))+
    xlab('Performance')+ ylab("\n ")+
    geom_errorbar(aes(ymin=lower, ymax=upper,col=Model),width=0.5,cex=1, position = position_dodge(0.6))+ 
    theme(plot.title=element_text(size=16,face="bold"),panel.background = element_blank(),
          panel.border = element_rect(fill=NA, colour="black"),
        axis.text.x=element_text(face="bold"),
        axis.title=element_text(size=12,face="bold"),
        strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold"),
        text=element_text(family="Century Gothic"))+
  scale_colour_manual(values=c("aquamarine3", "seagreen","salmon", "purple", "cornflowerblue", "darkgoldenrod2", "hotpink")) + scale_y_continuous(breaks = c(0,10,20,30,40,50,60,70,80,90,100)) +
    coord_flip()
p


p2<-ggplot(aes(x=Metric, y=Value, group=Model, colour=Model), data=Models_Performance) + geom_point(size=4, aes(color=Model), position=position_dodge(width=0.75)) + geom_segment(aes(xend=Metric, yend=0), position=position_dodge(width=0.75), size=1) +
  scale_colour_manual(values=c("aquamarine3", "seagreen","salmon", "purple", "cornflowerblue", "darkgoldenrod2", "hotpink")) + scale_y_continuous(breaks = c(0,10,20,30,40,50,60,70,80,90,100)) +  theme(plot.title=element_text(size=16,face="bold"),panel.background = element_blank(),panel.border = element_rect(fill=NA, colour="black"),
        axis.text.x=element_text(face="bold"),
        axis.title=element_text(size=12,face="bold"),
        strip.text.y = element_text(hjust=0,vjust = 1,angle=180,face="bold"),
        text=element_text(family="Segoe UI Historic", size=12))+ annotate("text", x=1, y=117, label="")+coord_flip()
ggsave(p2, dpi=400, height=7, width=8, filename = "C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Figure1.png")

kruskal.test(Models_Performance$Value[Models_Performance$Metric=="Accuracy"] ~ Models_Performance$Model[Models_Performance$Metric=="Accuracy"])

summary(aov(Models_Performance$Value[Models_Performance$Metric=="Sensitivity"] ~ Models_Performance$Model[Models_Performance$Metric=="Sensitivity"]))

```
    
    
    # P r o b a b i l i t y  P l o t
```{r}

modelProbs<-cbind(glmProbs[,1], rfProbs[,1], svmProbs[,1], nbProbs[,1], knnProbs[,1], annProbs[,1], glmnetProbs[,1], test.data12_i$dp12)
colnames(modelProbs)<-c("Logistic", "RF", "SVM", "NB", "KNN", "ANN", "GLMnet", "Outcome" )
modelProbs<-as.data.frame(modelProbs)
modelProbs$Outcome<-recode_factor(modelProbs$Outcome, "1"="Yes", "2"="No")
modelProbs

modelProbs<-melt(modelProbs, id.variables=c("Outcome"))

ggplot(modelProbs, aes(x=variable, y=value, fill=Outcome)) + geom_boxplot() + theme(panel.background = element_blank(),panel.border = element_rect(fill=NA, colour="black"),text=element_text(family="Segoe UI Historic", size=16)) + labs(y="Probability of Death/Progression\n", x="") + scale_y_continuous(labels=scales::percent, breaks = c(0, 0.10,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1)) + stat_compare_means(aes(label = sprintf("p = %5.4f", as.numeric(..p.format..))),method="wilcox.test") + scale_fill_manual(values=c("darkgoldenrod2", "azure4")) + annotate("text", x=1, y=1.1, label="", family="Segoe UI Historic") 




summary(modelProbs$value[modelProbs$Outcome=="Yes" & modelProbs$variable=="RF"])
summary(modelProbs$value[modelProbs$Outcome=="No" & modelProbs$variable=="RF"])
```
    

    
    # Descriptive Tables
```{r}
#Table1 -----------------------------------------------------------------------------------
load("C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/02-Data/20-05-2020/clinical_analysis_enrollment_2020-05-20.Rda")
load("C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/02-Data/20-05-2020/clab_analysis_2020-05-20.Rda")

data<-merge(clinical_analysis_enrollment[, c("record_id","redcap_data_access_group", "age_at_recruitment", "birth_sex", "weightkg","length","WAZ", "age_at_diagnosis", "age_at_art", "children_baseline_ART")], clab_analysis[which(clab_analysis$redcap_event_name=="enrollment"), c("record_id", "vl_rna_0", "cd4_p_0", "mother_SLEHI", "motart_adherence", "mother_cd4_count", "mother_vl_count", "death", "progression","dp6","dp12","dp", "dp_time", "t_fup")], by="record_id", all=T)

data<-data[data$record_id %in% mydata_dp12$record_id, ]
data$death<-factor(data$death)
data$progression<-factor(data$progression)
data$dp12<-mydata_dp12$dp12
summary(data)
data$WAZ<-as.numeric(data$WAZ$waz)
data$age_at_diagnosis<-data$age_at_diagnosis/30.5
data$age_at_art<-data$age_at_art/30.5
data$mother_SLEHI<-factor(data$mother_SLEHI)

tab<-compareGroups(dp12~., data=data, method=4)
table1<-createTable(tab, show.ratio = F, show.all = T)
export2word(table1, file="C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/Table1_20-05-2020.docx")

# --------------------------------------------------------------------------------------


#Table 2, testing and training
train.data12_i$set<-"training set"
test.data12_i$set<-"testing set"
data_total<-rbind(train.data12_i, test.data12_i)
tab<-compareGroups(set~., data=data_total, method=4)
table1<-createTable(tab, show.ratio = F, show.all = T)
export2word(table1, file="C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Table2_26-05-2020.docx")


# --------------------------------------------------------------------------------------

#imputed vs not imputed
train.data12$imputed<-"No"
colnames(train.data12_i)[1]<-"record_id"
train.data12_i$imputed<-"Yes"
mydata<-rbind(train.data12, train.data12_i)

tab<-compareGroups(imputed ~., data=mydata, method=4)
table1<-createTable(tab, show.ratio = F, show.all = T)
export2word(table1, file="C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Tables/SupTable1_train.docx")


#imputed vs not imputed
test.data12$imputed<-"No"
colnames(test.data12_i)[1]<-"record_id"
test.data12_i$imputed<-"Yes"
mydata_test<-rbind(test.data12, test.data12_i)

tab<-compareGroups(imputed ~., data=mydata_test, method=4)
table1<-createTable(tab, show.ratio = F, show.all = T)
export2word(table1, file="C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Tables/SupTable1_test.docx")


```
    
    
   
   
    
    ##   C O M P L E T E   C A S E S
    
     
```{r}
summary(mydata_dp12_completecase)
str(mydata_dp12_completecase)

mydata_dp12_completecase$birth_sex<-recode_factor(mydata_dp12_completecase$birth_sex, "1"="Female", "2"="Male")
mydata_dp12_completecase$dp12<-recode_factor(mydata_dp12_completecase$dp12, "1"="Yes", "0"="No")

```
     
 
```{r}

#train and test
#training
train.rows<- createDataPartition(y= mydata_dp12_completecase$dp12, p=0.7, list = FALSE)
train.data12<- mydata_dp12_completecase[train.rows,] # 70% data goes in here
summary(factor(train.data12$dp12))

#testing
test.data12<- mydata_dp12_completecase[-train.rows,] # 30% data goes in here
summary(factor(test.data12$dp12))


save(train.data12,file = paste(paste("train.data_dp12_cc",Sys.Date(),sep = "_"),".Rda",sep=""))
save(test.data12,file = paste(paste("test.data_dp12_cc",Sys.Date(),sep = "_"),".Rda",sep=""))

```
    

    # Model 1: Logistic Regression

```{r}
mod<-glm(dp12 ~ .,data = train.data12_completecase[,c(2:13)], family=binomial)
stepAIC(mod)


glmFit <- train(dp12 ~ .,data = train.data12_completecase[,c(2:13)],method = "glm", family="binomial",trControl=ctrl)

glmClasses <- predict(glmFit, newdata = test.data12_completecase)
glmProbs <- predict(glmFit, newdata = test.data12_completecase, type = "prob")

confusionMatrix(data = glmClasses , test.data12_completecase$dp12) #sensitiviy 0.67 #specificity 0.81 

glmFit_cc<-glmFit
saveRDS(glmFit_cc, "C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Models/completecases/glmFit_cc.rds")

```


        # Model 2: Random Forest
```{r}
rfFitcc <- train(dp12 ~ .,data = train.data12_completecase[,c(2:13)],method = "rf",trControl=ctrl, tuneLength=10)
rdaGrid = data.frame(mtry = 3)
rfFitcc <- train(dp12 ~ .,data = train.data12_completecase[,c(2:13)],method = "rf",trControl=ctrl, tuneGrid = rdaGrid)

rfClasses <- predict(rfFitcc, newdata = test.data12_completecase)
rfProbs <- predict(rfFitcc, newdata = test.data12_completecase, type = "prob")

confusionMatrix(data = rfClasses , test.data12_completecase$dp12) #sensitiviy 0.7 #specificity 0.8 

saveRDS(rfFitcc, "C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Models/completecases/rfFitcc.rds")

```
        

        # Model 3: Support Vector Machine
```{r}
svmFitcc <- train(dp12 ~ .,data = train.data12_completecase[,c(2:13)],method = "svmRadial",trControl=ctrl, tuneLength=10)
rdaGrid = data.frame(C = 4, sigma=1.899027e-10)
svmFitcc <- train(dp12 ~ .,data = train.data12_completecase[,c(2:13)],method = "svmRadial",trControl=ctrl, tuneGrid = rdaGrid)

svmClasses <- predict(svmFitcc, newdata = test.data12_completecase)
svmProbs <- predict(svmFitcc, newdata = test.data12_completecase, type = "prob")

confusionMatrix(data = svmClasses , test.data12_completecase$dp12) #sensitiviy 0.8 #specificity 0.7

saveRDS(svmFitcc, "C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Models/completecases/svmFitcc.rds")

```
        

          # Model 4: Naive Bayes
```{r}
nbFitcc <- train(dp12 ~ .,data = train.data12_completecase[,c(2:7,9:13)],method = "nb")
rdaGrid = data.frame(fL=0, usekernel=TRUE, adjust=1)
nbFitcc <- train(dp12 ~ .,data = train.data12_completecase[,c(2:7,9:13)],method = "nb",trControl=ctrl, tuneGrid = rdaGrid)

nbClasses <- predict(nbFitcc, newdata = test.data12_completecase)
nbProbs <- predict(nbFitcc, newdata = test.data12_completecase, type = "prob")

confusionMatrix(data = nbClasses , test.data12_completecase$dp12) #sensitiviy 0.8 #specificity 0.7

saveRDS(nbFitcc, "C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Models/completecases/nbFitcc.rds")
```
        
        # Model 5: KNN
```{r}
knnFitcc <- train(dp12 ~ .,data = train.data12_completecase[,c(2:13)],method = "knn",trControl=ctrl, tuneLength=10)
rdaGrid = data.frame(k=5)
knnFit <- train(dp12 ~ .,data = train.data12_completecase[,c(2:13)],method = "knn",trControl=ctrl, tuneGrid = rdaGrid)

knnClasses <- predict(knnFitcc, newdata = test.data12_completecase)
knnProbs <- predict(knnFitcc, newdata = test.data12_completecase, type = "prob")

confusionMatrix(data = knnClasses , test.data12_completecase$dp12) #sensitiviy 0.8 #specificity 0.7

saveRDS(knnFitcc , "C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Models/completecases/knnFitcc .rds")
```
        
    # Model 6: ANN
```{r}
annFitcc <- train(dp12 ~ .,data = train.data12_completecase[,c(2:13)],method = "nnet",trControl=ctrl, tuneLength=10)
rdaGrid = data.frame(size=15, decay=0.042)
annFitcc <- train(dp12 ~ .,data = train.data12_completecase[,c(2:13)],method = "nnet",trControl=ctrl, tuneGrid = rdaGrid)

annClasses <- predict(annFitcc, newdata = test.data12_completecase)
annProbs <- predict(annFitcc, newdata = test.data12_completecase, type = "prob")

confusionMatrix(data = annClasses , test.data12_completecase$dp12) #sensitiviy 0.8 #specificity 0.7

saveRDS(annFitcc , "C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Models/completecases/annFitcc.rds")
```
   
    # Model 7: GLMNET
```{r}
glmnetFitcc <- train(dp12 ~ .,data = train.data12_completecase[,c(2:13)],method = "glmnet",trControl=ctrl, tuneLength=10)
rdaGrid = data.frame(alpha=0.5, lambda=0.001)
glmnetFitcc <- train(dp12 ~ .,data = train.data12_completecase[,c(2:13)],method = "glmnet",trControl=ctrl, tuneGrid = rdaGrid)

glmnetClasses <- predict(glmnetFitcc , newdata = test.data12_completecase)
glmnetProbs <- predict(glmnetFitcc , newdata = test.data12_completecase, type = "prob")

confusionMatrix(data = glmnetClasses , test.data12_completecase$dp12) #sensitiviy 0.8 #specificity 0.7

saveRDS(glmnetFitcc, "C:/Users/sarad/OneDrive/Escritorio/EPIICAL/EARTH/19- Machine Learning/ML May 2020/Models/completecases/glmnetFitcc.rds")


```