---
title: "Big Data Machine Learning Methods for Metabolomics"
subtitle: Part II
author: 
  - Kangni Alemdjrodo, Department of Statistics, Purdue University
  - Dabao Zhang, Department of Statistics, Purdue University
  - Daniel Raftery, Department of Anesthesiology and Pain Medicine, University of Washington
  - Min Zhang, Department of Statistics, Purdue University
date: "Summer 2022"
output:
  pdf_document:
    number_sections: yes
    toc: yes
    toc_depth: 3
  html_document:
    number_sections: yes
    toc: yes
    toc_depth: 3
---


# Introduction

This is the second part of the notebook. We have to reload the packages and data into R.

```{r, message=FALSE, warning=FALSE}
library(caret)   ## For model tuning using resampling
library(dplyr)   ## For some data preparation
library(Matrix)  ## For sparse matrix
library(missRanger) ## For missing value imputation
library(mixOmics)## For pls-da
library(POCRE)   ## For gps
library(ROCR)    ## For ROC curves
library(SIS)     ## For sis
library(splitTools) ## For stratified samples
library(xgboost) ## For xgboost

# Load the data sets
targeted  <- read.csv("targeted.csv") 
subject_info <- read.csv("subject_info.csv")

# Reclassify "CM", "CO", "CY" as "Control"
subject_info[subject_info$Type=="CM",]$Type <- "Control"
subject_info[subject_info$Type=="CO",]$Type <- "Control"
subject_info[subject_info$Type=="CY",]$Type <- "Control"

# Exclude PD patients
targeted <- targeted[subject_info$Type!="PD",]
subject_info <- subject_info[subject_info$Type!="PD",]

# Type: 0-Contro; 1-AD
subject_info$Type <- factor(subject_info$Type,levels=c("Control","AD"),labels=c(0,1))
# Sex: 0-M; 1-F
subject_info$Sex <- factor(subject_info$Sex,levels=c("M","F"),labels=c(0,1))

# Impute the missing values
set.seed(20220620) ## To replicate the results
targeted <- missRanger(targeted[, -1], num.trees = 100) ## Remove id before imputation
```

You may practice what you have learned in the previous part of notebook by building predictive models to predict AD patients using both demographic and metabolomic information. For this purpose, we can first create `csf_md`, which include both disease type (control vs. AD), demographic information (*Age* and *Sex*), and metabolites levels. 

```{r}
csf_md <- data.frame(Type=subject_info$Type, targeted, Age=subject_info$Age, Sex=subject_info$Sex)
```

You can replace `csf_m` by `csf_md` in the previous part of the notebook, and then rerun all the codes.



# Predictive Models Including Interactions between Metabolomic and Demographic Varaibles 

Here we intend to build predictive models to predict AD patients using both demographic and metabolomic information as well as their interactions. We first create `csf_mxd`, which include both disease type (control vs. AD), demographic inforamtion, metabolites levels, and the interactions between demographic and metabolomic variables. 


```{r}
## Standardize both metabolites and demographics variables
metab <- scale(targeted)
subject_info$Sex <- as.numeric(subject_info$Sex)
subject_info$Age <- as.numeric(subject_info$Age)
demog <- scale(subject_info[,3:4])

##  Create design matrix including metabolite*demographic variable interactions
design <- model.matrix(Type~metab*demog-1, data=subject_info)

## Create the corresponding data frame
csf_mxd <- data.frame(Type=subject_info$Type, design)

dim(csf_mxd)
```

We created a total of 216 metabolites-demographics interaction terms. We can repeat the analysis on `csf_m` by replacing `csf_md by `csf_mxd`.



###  Logistic regression

```{r, warning=FALSE}
set.seed(20220620) ## To replicate the results
start.time <- Sys.time() # Start the clock!
nsim <- 100 ## Number of repetitions

csf_mxd.logit.test.preds <- vector(mode="list", length=nsim)
csf_mxd.logit.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  train.idx <- partition(csf_mxd$Type, p=c(train=0.8, test=0.2)) ## library(SpliTools)
  train <- csf_mxd[train.idx$train,]
  test  <- csf_mxd[train.idx$test,]
  
  ## Fit LOGISTIC 
  csf_mxd.logit.train <- glm(Type~., data=train, family=binomial(link="logit"))

  ## Predict
  csf_mxd.logit.test.preds[[j]] <- predict(csf_mxd.logit.train, test, type="response")
  csf_mxd.logit.test.labels[[j]]<- as.numeric(as.character(test$Type))
}

csf_mxd.logit.ROCR <- ROCR::prediction(csf_mxd.logit.test.preds, csf_mxd.logit.test.labels)
csf_mxd.logit.roc <- performance(csf_mxd.logit.ROCR,"tpr","fpr")
csf_mxd.logit.auc <- performance(csf_mxd.logit.ROCR, measure = "auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```


Let us plot the ROC curve and get the AUC.

```{r}
mar.default <- c(5,4,4,2) + 0.1
par(mar = mar.default + c(0, 2, 0, 0)) 
plot(csf_mxd.logit.roc,col="grey50",lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_mxd.logit.roc,lwd=3,avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="Logistic Regression with All Predictors")
text(x=0.8, y=0.1,paste("AUC=", round(median(unlist(csf_mxd.logit.auc@y.values)),2)),cex = 1.5)  # AUC
abline(a=0, b=1) 
par(new=FALSE)
```

While the above plot shows with AUC=0.67, we can further confirm this and calculate its standard error.

```{r}
round(median(unlist(csf_mxd.logit.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_mxd.logit.auc@y.values)),2)/10  ## standard error
```



### PLS-DA 

A PLS-DA model is fitted with 10 components to evaluate the performance and the number of components necessary for the final model (see below).

```{r}
csf_mxd.plsda <- plsda(csf_mxd[,-1], csf_mxd$Type, ncomp=10)  # set ncomp to 10
```

The samples are then projected into the subspace spanned by the first two components.

```{r}
plotIndiv(csf_mxd.plsda, comp=1:2, group=csf_mxd$Type, ind.names=FALSE,
          ellipse=TRUE, legend=TRUE, title='PLS-DA with All Predictors')
```

From the sample plot, we observe a better separation of 2 groups compared to an unsupervised PCA sample plot (skipped). Confidence ellipses for each class are plotted to highlight the strength of the discrimination (confidence level set to 95% by default, see argument ellipse.level). 

The prediction area can be visualized by calculating a background surface first, before overlaying the sample plot. See `? background.predict` for more details. The prediction distance needs to be specified:

```{r}
bg.csf_mxd.plsda <- background.predict(csf_mxd.plsda, comp.predicted=2, dist="max.dist")
plotIndiv(csf_mxd.plsda, comp=1:2, group=csf_mxd$Type, ind.names=FALSE, title="Maximum Distance",
          legend=TRUE,  background=bg.csf_mxd.plsda)
```


**Performance**

The classification performance of the PLS-DA model is assessed with the `perf` function using 5-fold
cross-validation repeated 50 times. The number of repeats is necessary to ensure a good estimate of the classification error rate (as the CV-folds are determined in a random manner). From the performance results we can decide on the choice of the number of components of the final PLSDA model. The maximal number of components assessed corresponds to the number of components of the `csf_mxd.plsda` model ran above.

This takes a couple of minutes (2 min) to run.

```{r}
set.seed(20220620) # for reproducibility, only when the `cpus' argument is not used
start.time <- Sys.time() # Start the clock!

csf_mxd.plsda.perf <- perf(csf_mxd.plsda, validation="Mfold", folds=5,
                         progressBar=FALSE, auc=TRUE, nrepeat=50)

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken) 
```

**Error rates**
```{r}
plot(csf_mxd.plsda.perf, col=color.mixo(5:7), sd=TRUE, legend.position="horizontal")
```


From the performance plot, we observe that the centroid error rate decreases slowly from 1 to 2 components. The error rates does not stabilize after 2 components

The perf function outputs the optimal number of components in `csf_mxd.plsda.perf$choice.ncomp` based on t-tests that test for a significant difference in the mean error rate between components. The result is output for each performance measure and distance, when applicable.

```{r}
ncomps.csf_mxd.plsda <- csf_mxd.plsda.perf$choice.ncomp
ncomps.csf_mxd.plsda
```

```{r}
csf_mxd.plsda.perf$error.rate
```

**Conclusion:** The performance displayed suggests that `ncomp = 1`  with centroids distance is sufficient to achieve performance (0.31 error rate).


#### PLS-DA final model performance

```{r}
ncomps.csf_mxd.plsda = 1 # Number of component for PLS-DA 
set.seed(20220620) ## To replicate the results
start.time <- Sys.time() # Start the clock!
nsim <- 100 ## Number of repetitions
csf_mxd.plsda.test.preds <- vector(mode="list", length=nsim)
csf_mxd.plsda.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  train.idx <- partition(csf_mxd$Type, p=c(train=0.8, test=0.2))
  train <- csf_mxd[train.idx$train,]
  test  <- csf_mxd[train.idx$test,]
  
  ## Fit PLS-DA
  csf_mxd.plsda.train <- plsda(train[,-1], train$Type, ncomp=ncomps.csf_mxd.plsda)
  
  ## Predict
  csf_mxd.plsda.test.preds[[j]] <- predict(csf_mxd.plsda.train, test[,-1], 
                                         dist="centroids.dist")$predict[,,ncomps.csf_mxd.plsda][,2]
  csf_mxd.plsda.test.labels[[j]]<- as.numeric(test$Type)
}

csf_mxd.plsda.ROCR = ROCR::prediction(csf_mxd.plsda.test.preds, csf_mxd.plsda.test.labels)
csf_mxd.plsda.roc = performance(csf_mxd.plsda.ROCR,"tpr","fpr")
csf_mxd.plsda.auc = performance(csf_mxd.plsda.ROCR, measure="auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```


**ROC curve and AUC**

```{r}
mar.default <- c(5,4,4,2) + 0.1
par(mar = mar.default + c(0, 2, 0, 0)) 
plot(csf_mxd.plsda.roc,col="grey50", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_mxd.plsda.roc, lwd=3, avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="PLS-DA with All Predictors")
text(x=0.8, y=0.1, paste("AUC=", round(median(unlist(csf_mxd.plsda.auc@y.values)),2)), cex=1.5)  # AUC
abline(a=0, b=1)
par(new=FALSE)
```

From the above figure, we see that PLS-DA could build predictive models with All Predictors with AUC at 0.72, which is confirmed in the following with standard error.

```{r}
round(median(unlist(csf_mxd.plsda.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_mxd.plsda.auc@y.values)),2)/10  ## standard error
```



### XGBoost

We will take XGBoost to build tree ensemble models (XGBoost-T) and linear-regression ensemble models (XGBoost-L), respectively.


#### XGBoost-T

##### Determine optimal parameters to use with the XGBoost-T

**XGBoost-T with default parameters**

```{r}
start.time <- Sys.time() # Start the clock!
set.seed(20220620)

# Split the data
train.idx <- partition(csf_mxd$Type, p=c(train=0.8, test=0.2))
train <- csf_mxd[train.idx$train,]
test  <- csf_mxd[train.idx$test,]
  
### Make sure the dependent variable is numeric 0,1,2... (required by XGboost package)
train.X <- sparse.model.matrix(Type ~ .-1, data = train)
train.Y <- as.numeric(as.character(train$Type))

test.X <- sparse.model.matrix(Type~.-1, data = test)
test.Y <- as.numeric(as.character(test$Type))
  
## Convert the train and test data into xgb matrix type (required by XGboost package)
xgb_train = xgb.DMatrix(data=train.X, label=train.Y)
xgb_test = xgb.DMatrix(data=test.X, label=test.Y)
  
## Defining the parameters
xgb_params <- list("objective"="binary:logistic", ## logistic regression
                   "eval_metric"="logloss") ## multiclass logloss (used in classification)
watchlist <- list(train=xgb_train, test=xgb_test)

xgbc = xgb.train(booster="gbtree",
                 params=xgb_params,
                 data=xgb_train,
                 max_depth = 6,
                 eta = 0.3,
                 gamma=0,
                 nthread = 1,
                 nrounds=100, ## number of XGBoost iterations
                 watchlist=watchlist,
                 min_child_weight=1, ##
                 subsample=1, 
                 colsample_bytree=1, ##
                 verbose=F) ## Print or not

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
# print(xgbc)
```

**Check error rates**

Visually

```{r}
# Training & test error plot
error_xgbc <- data.frame(xgbc$evaluation_log)
plot(error_xgbc$iter, error_xgbc$train_logloss, col='blue', type="l", ylim=c(0,1))
lines(error_xgbc$iter, error_xgbc$test_logloss, col='red')
legend(x=60,y=0.4,c("Train","Test"), col=c('blue','red'), lty=1)
```


Using the Root mean square error(RMSE) 

```{r}
csf_mxd.xgbtree.test.preds <- as.numeric(predict(xgbc, xgb_test))
csf_mxd.xgbtree.test.labels <- test.Y
rmse.xgbtree <- sqrt(mean((csf_mxd.xgbtree.test.preds-csf_mxd.xgbtree.test.labels)^2)) 
rmse.xgbtree # RMSE
```

**Tune parameters**

This tunning process can be time consuming (here 2 min)

```{r, message=FALSE, warning=FALSE}
set.seed(20220620)
start.time <- Sys.time() # Start the clock!

xgbtree_caret <- train(x=as.matrix(train.X) ,
                   y=as.factor(train.Y),
                   method="xgbTree",
                   trControl=trainControl(method="repeatedcv",
                                          number=3,
                                          repeats=2,
                                          verboseIter=FALSE),
                   tuneGrid=expand.grid(nrounds=c(50,100,500),
                                        eta=c(0.01,0.05,0.1),
                                        max_depth=6,
                                        colsample_bytree=c(0.5,1),
                                        subsample=c(0.5,1),
                                        gamma=0,
                                        min_child_weight=c(0.5,1,2)))

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```

The best parameters

```{r}
xgbtree_caret$bestTune 
```


Run the XGBoost with the tuned parameters

```{r}
 set.seed(20220620)
 xgbtree_tuned <- xgb.train(booster="gbtree",
                    params=xgb_params,
                    data=xgb_train,
                    max_depth=6,
                    eta=xgbtree_caret$bestTune$eta,
                    gamma=0,
                    nthread=1,
                    nrounds=xgbtree_caret$bestTune$nrounds, ## number of XGBoost iterations
                    watchlist=watchlist,
                    min_child_weight=xgbtree_caret$bestTune$min_child_weight, ##
                    subsample=xgbtree_caret$bestTune$subsample, 
                    colsample_bytree=xgbtree_caret$bestTune$colsample_bytree, ##
                    verbose=F) ## Print or not
print(xgbtree_tuned)
```

**Check error rates**

Visually

```{r}
# Training & test error plot
error_xgbtree_tuned <- data.frame(xgbtree_tuned$evaluation_log)
plot(error_xgbtree_tuned$iter, error_xgbtree_tuned$train_logloss, col='blue', type="l", ylim=c(0,1))
lines(error_xgbtree_tuned$iter, error_xgbtree_tuned$test_logloss, col='red')
legend(x=40,y=0.9,c("Train","Test"), col=c('blue','red'), lty=1)
```

Using the Root mean square error(RMSE) 

```{r}
csf_mxd.xgbtree.test.preds <- as.numeric(predict(xgbtree_tuned, xgb_test))
csf_mxd.xgbtree.test.labels <- test.Y
rmse.xgbtree_tuned <- sqrt(mean((csf_mxd.xgbtree.test.preds-csf_mxd.xgbtree.test.labels)^2))
rmse.xgbtree_tuned # RMSE
```


**Note**: The RMSE was 0.489 with the default parameters.



##### Final model performance

```{r}
set.seed(20220620) ## To replicate the results
nsim <- 100

csf_mxd.xgbtree.test.preds <- vector(mode="list", length=nsim)
csf_mxd.xgbtree.test.labels <- vector(mode="list", length=nsim)

## Extract x - feature and y - label parts.
## To use the XGBoost algorithm we need to create a matrix and use the one-hot encoding.
## This will create a dummy variable for the factor variables: DONE

for(j in 1:nsim) {
  train.idx <- partition(csf_mxd$Type, p=c(train=0.8, test=0.2))
  
  train <- csf_mxd[train.idx$train,]
  test  <- csf_mxd[train.idx$test,]
  
  ### Make sure the dependent variable is numeric 0,1,2... (required by XGboost package)
  train.X <- sparse.model.matrix(Type~.-1, data=train)
  train.Y <- as.numeric(as.character(train$Type))
  
  test.X <- sparse.model.matrix(Type~.-1, data=test)
  test.Y <- as.numeric(as.character(test$Type))
  
  ## Convert the train and test data into xgb matrix type (required by XGboost package)
  xgb_train = xgb.DMatrix(data=train.X, label=train.Y)
  xgb_test = xgb.DMatrix(data=test.X, label=test.Y)
  
  ## Defining the parameters
  xgb_params <- list("objective"="binary:logistic", ## logistic regression for binary classification
                     "eval_metric" = "logloss") ## multiclass logloss (used in classification)
  watchlist <- list(train=xgb_train, test=xgb_test)
  
  xgbc <- xgb.train(booster="gbtree",
                    params=xgb_params,
                    data=xgb_train,
                    max_depth=6,
                    eta=xgbtree_caret$bestTune$eta,
                    gamma=0,
                    nthread=1,
                    nrounds=xgbtree_caret$bestTune$nrounds, ## number of XGBoost iterations
                    watchlist=watchlist,
                    min_child_weight=xgbtree_caret$bestTune$min_child_weight, ##
                    subsample=xgbtree_caret$bestTune$subsample, 
                    colsample_bytree=xgbtree_caret$bestTune$colsample_bytree, ##
                    verbose=F) ## Print or not
  ### print(xgbc)

  ## The model is ready and we can predict the test data.
  ## Predicting test data (Hold-out)
  csf_mxd.xgbtree.test.preds[[j]] <- as.numeric(predict(xgbc, xgb_test))
  csf_mxd.xgbtree.test.labels[[j]] <- test.Y
}

csf_mxd.xgbtree.ROCR = ROCR::prediction(csf_mxd.xgbtree.test.preds, csf_mxd.xgbtree.test.labels)
csf_mxd.xgbtree.roc = performance(csf_mxd.xgbtree.ROCR,"tpr","fpr")
csf_mxd.xgbtree.auc = performance(csf_mxd.xgbtree.ROCR, measure="auc")
```


ROC curve and AUC

```{r}
mar.default <- c(5,4,4,2) + 0.1   # plot margins
par(mar = mar.default + c(0, 2, 0, 0)) 

plot(csf_mxd.xgbtree.roc,col="grey50", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_mxd.xgbtree.roc, lwd=3, avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="XGBoost-T with All Predictors")
text(x=0.8, y=0.1, paste("AUC=", round(median(unlist(csf_mxd.xgbtree.auc@y.values)),2)), cex=1.5) # AUC
abline(a=0, b=1)
par(new=FALSE)
```


```{r}
round(median(unlist(csf_mxd.xgbtree.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_mxd.xgbtree.auc@y.values)),2)/10  ## standard error
```


####  XGBoost-L
##### Determine optimal parameters to use with the XGBoost-L


**XGBoost-L with default parameters**


```{r}
set.seed(20220620)
start.time <- Sys.time() # Start the clock!
train.idx <- partition(csf_mxd$Type, p=c(train=0.8, test=0.2)) ## library(SpliTools)
  
train <- csf_mxd[train.idx$train,]
test  <- csf_mxd[train.idx$test,]

### Make sure the dependent variable is numeric 0,1,2... (required by XGboost package)
train.X <- as.matrix(train[,-1])
train.Y <- as.numeric(as.character(train$Type))
  
test.X <- as.matrix(test[,-1])
test.Y <- as.numeric(as.character(test$Type))
  
## Need to convert the train and test data into xgb matrix type (required by XGboost package)
xgb_train <- xgb.DMatrix(data=train.X, label=train.Y)
xgb_test <- xgb.DMatrix(data=test.X, label=test.Y)

## Defining the parameters
#numberOfClasses <- length(unique(csf_mxd$Type))
xgb_params <- list("objective"="binary:logistic",  ## logistic regression for binary classification
                   "eval_metric"="logloss") ## multiclass logloss (used in classification)
watchlist <- list(train=xgb_train, test=xgb_test)
  
xgbc <- xgb.train(data=xgb_train, booster="gblinear",
                  params=xgb_params,
                  nrounds=100,
                  verbose=0,
                  watchlist=watchlist)

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken) 
```


**Check error rates**

```{r}
csf_mxd.xgblin.test.preds <- as.numeric(predict(xgbc, xgb_test))
csf_mxd.xgblin.test.labels <- test.Y
rmse_xgbc <- sqrt(mean((csf_mxd.xgblin.test.preds-csf_mxd.xgblin.test.labels)^2))
rmse_xgbc # RMSE .
```


**Tune parameters**

```{r, message=FALSE, warning=FALSE}
set.seed(20220620)
start.time <- Sys.time() # Start the clock!

xgblin_caret <- train(x=as.matrix(train.X) ,
                   y=as.factor(train.Y),
                   method="xgbLinear",
                   trControl=trainControl(method="repeatedcv",
                                          number=5,
                                          repeats=10,
                                          verboseIter=FALSE),
                   tuneGrid=expand.grid(nrounds=c(1,2,5),
                                        lambda =0,
                                        alpha = 1,
                                        eta=0.3))

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```

The best parameters are
```{r}
xgblin_caret$bestTune 
```


```{r}
xgblin_tuned <- xgb.train(data=xgb_train, booster="gblinear",
                        params=xgb_params,
                        nrounds=xgblin_caret$bestTune$nrounds,
                        verbose=0,
                        watchlist=watchlist)
```

```{r}
csf_mxd.xgblin.test.preds_tuned <- as.numeric(predict(xgblin_tuned, xgb_test))
csf_mxd.xgblin.test.labels_tuned <- test.Y
rmse_xgblin_tuned <- sqrt(mean((csf_mxd.xgblin.test.preds_tuned-csf_mxd.xgblin.test.labels_tuned)^2)) 
rmse_xgblin_tuned # RMSE
```

The error rate with default parameters was 0.586.


#####  Final model performance

```{r}
set.seed(20220620) ## To replicate the results
start.time <- Sys.time() # Start the clock!
nsim <- 100

csf_mxd.xgblin.test.preds <- vector(mode="list", length=nsim)
csf_mxd.xgblin.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  ## Extract x - feature and y - label parts.
  ## To use the XGBoost algorithm we need to create a matrix and use the one-hot encoding.
  ## This will create a dummy variable for the factor variables: DONE
  
  train.idx <- partition(csf_mxd$Type, p=c(train=0.8, test=0.2)) ## library(SpliTools)
  
  train <- csf_mxd[train.idx$train,]
  test  <- csf_mxd[train.idx$test,]
  
  train.X <- as.matrix(train[,-1])
  train.Y <- as.numeric(as.character(train$Type))
  
  test.X <- as.matrix(test[,-1])
  test.Y <- as.numeric(as.character(test$Type))
  
  ## Convert the train and test data into xgb matrix type (required by XGboost package)
  xgb_train <- xgb.DMatrix(data=train.X, label=train.Y)
  xgb_test <- xgb.DMatrix(data=test.X, label=test.Y)

  ## Defining the parameters
  xgb_params <- list("objective"="binary:logistic",  ## logistic regression
                     "eval_metric"="logloss") ## multiclass logloss (used in classification)
  watchlist <- list(train=xgb_train, test=xgb_test)
  
  xgbc <- xgb.train(data=xgb_train, booster="gblinear",
                    params=xgb_params,
                    nrounds=xgblin_caret$bestTune$nrounds,
                    verbose=0,
                    watchlist=watchlist)

  ## The model is ready and we can predict the test data.
  ## Predicting test data (Hold-out)
  csf_mxd.xgblin.test.preds[[j]] <- as.numeric(predict(xgbc, xgb_test))
  csf_mxd.xgblin.test.labels[[j]]<- as.numeric(test.Y)
}

csf_mxd.xgblin.ROCR = ROCR::prediction(csf_mxd.xgblin.test.preds, csf_mxd.xgblin.test.labels)

csf_mxd.xgblin.roc = performance(csf_mxd.xgblin.ROCR,"tpr","fpr")
csf_mxd.xgblin.auc = performance(csf_mxd.xgblin.ROCR, measure="auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken) 
```



**ROC curve and AUC**

```{r}
mar.default <- c(5,4,4,2) + 0.1
par(mar = mar.default + c(0, 2, 0, 0)) 
plot(csf_mxd.xgblin.roc, col="grey50", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_mxd.xgblin.roc,lwd=3,avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="XGBoost-L with All Predictors")
abline(a=0, b=1) 
text(x=0.8, y=0.1,paste("AUC=", round(median(unlist(csf_mxd.xgblin.auc@y.values)),2)),cex = 1.5)  # Show AUC
par(new=FALSE)
```


```{r}
round(median(unlist(csf_mxd.xgblin.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_mxd.xgblin.auc@y.values)),2)/10  ## standard error
```



##  Predictive Models with Selected Predictors

We can consider two different methods to select a pre-specified number of predictors to build the predictive models for AD patients. GPS & SIS...

Note that a cross-validation method can be used to optimize the number of metabolites to select.

### Logistic regression with 35 Predictors Selected by GPS


```{r, warning=FALSE}
set.seed(20220620) ## To replicate the results
start.time <- Sys.time() # Start the clock!
nsim <- 100 ## Number of repetitions

csf_mxd35_GPS.logit.test.preds <- vector(mode="list", length=nsim)
csf_mxd35_GPS.logit.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  train.idx <- partition(csf_mxd$Type, p = c(train=0.8, test=0.2)) ## library(SpliTools)
  train <- csf_mxd[train.idx$train,]
  test  <- csf_mxd[train.idx$test,]
  
  m35_GPS <- gps(train.Y, train.X, family='binomial', maxvar=35, maxcmp=5)

  ## Fit LOGISTIC 
  csf_mxd35_GPS.logit.train <- glm(Type~., data=train[,m35_GPS], family=binomial(link="logit"))

  ## Predict
  csf_mxd35_GPS.logit.test.preds[[j]] <- predict(csf_mxd35_GPS.logit.train, test[,m35_GPS], type="response")
  csf_mxd35_GPS.logit.test.labels[[j]]<- as.numeric(as.character(test$Type))
}

csf_mxd35_GPS.logit.ROCR <- ROCR::prediction(csf_mxd35_GPS.logit.test.preds, csf_mxd35_GPS.logit.test.labels)
csf_mxd35_GPS.logit.roc <- performance(csf_mxd35_GPS.logit.ROCR,"tpr","fpr")
csf_mxd35_GPS.logit.auc <- performance(csf_mxd35_GPS.logit.ROCR, measure = "auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```


Let us plot the ROC curve and get the AUC.

```{r}
mar.default <- c(5,4,4,2) + 0.1
par(mar = mar.default + c(0, 2, 0, 0)) 
plot(csf_mxd35_GPS.logit.roc,col="grey50",lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_mxd35_GPS.logit.roc,lwd=3,avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="Logistic Regression with 35 Predictors by GPS")
text(x=0.8, y=0.1,paste("AUC=", round(median(unlist(csf_mxd35_GPS.logit.auc@y.values)),2)),cex = 1.5)  # AUC
abline(a=0, b=1) 
par(new=FALSE)
```

While the above plot shows with AUC=0.67, we can further confirm this and calculate its standard error.

```{r}
round(median(unlist(csf_mxd35_GPS.logit.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_mxd35_GPS.logit.auc@y.values)),2)/10  ## standard error
```


### Logistic regression with 35 Predictors Selected by SIS


```{r, warning=FALSE}
set.seed(20220620) ## To replicate the results
start.time <- Sys.time() # Start the clock!
nsim <- 100 ## Number of repetitions

csf_mxd35_SIS.logit.test.preds <- vector(mode="list", length=nsim)
csf_mxd35_SIS.logit.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  train.idx <- partition(csf_mxd$Type, p = c(train=0.8, test=0.2)) ## library(SpliTools)
  train <- csf_mxd[train.idx$train,]
  test  <- csf_mxd[train.idx$test,]
  
  m35_SIS <- SIS(train.X, train.Y, family='binomial', iter=FALSE, tune='cv', nsis=35)$ix0

  ## Fit LOGISTIC 
  csf_mxd35_SIS.logit.train <- glm(Type~., data=train[,m35_SIS], family=binomial(link="logit"))

  ## Predict
  csf_mxd35_SIS.logit.test.preds[[j]] <- predict(csf_mxd35_SIS.logit.train, test[,m35_SIS], type="response")
  csf_mxd35_SIS.logit.test.labels[[j]]<- as.numeric(as.character(test$Type))
}

csf_mxd35_SIS.logit.ROCR <- ROCR::prediction(csf_mxd35_SIS.logit.test.preds, csf_mxd35_SIS.logit.test.labels)
csf_mxd35_SIS.logit.roc <- performance(csf_mxd35_SIS.logit.ROCR,"tpr","fpr")
csf_mxd35_SIS.logit.auc <- performance(csf_mxd35_SIS.logit.ROCR, measure = "auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)

```


Let us plot the ROC curve and get the AUC.

```{r}
mar.default <- c(5,4,4,2) + 0.1
par(mar = mar.default + c(0, 2, 0, 0)) 
plot(csf_mxd35_SIS.logit.roc,col="grey50",lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_mxd35_SIS.logit.roc,lwd=3,avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="Logistic Regression with 35 Predictors by SIS")
text(x=0.8, y=0.1,paste("AUC=", round(median(unlist(csf_mxd35_SIS.logit.auc@y.values)),2)),cex = 1.5)  # AUC
abline(a=0, b=1) 
par(new=FALSE)
```

While the above plot shows with AUC=0.67, we can further confirm this and calculate its standard error.

```{r}
round(median(unlist(csf_mxd35_SIS.logit.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_mxd35_SIS.logit.auc@y.values)),2)/10  ## standard error
```



###  PLS-DA with 35 Predictors Selected by GPS

**Tune the number of components for the PLS-DA**

```{r}
set.seed(20220620) # for reproducibility, only when the `cpus' argument is not used
start.time <- Sys.time() # Start the clock!

train.idx <- partition(csf_mxd$Type, p = c(train=0.8, test=0.2)) ## library(SpliTools)
  
train <- csf_mxd[train.idx$train,]
test  <- csf_mxd[train.idx$test,]
  
train.X <- as.matrix(train[,-1])
train.Y<- as.numeric(as.character(train$Type))
  
m35_GPS <- gps(train.Y, train.X, family='binomial', maxvar=35, maxcmp=4)
m35_GPS.plsda <- plsda(train.X[, m35_GPS], train.Y, ncomp=10)  # set ncomp to 10


m35_GPS.plsda.perf <- perf(m35_GPS.plsda, validation="Mfold", folds=5,
                         progressBar=FALSE, auc=TRUE, nrepeat=50)
                         
ncomps.m35_GPS.plsda <- m35_GPS.plsda.perf$choice.ncomp
ncomps.m35_GPS.plsda

m35_GPS.plsda.perf$error.rate

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken) 
```

**Final model**

```{r, message=FALSE, warning=FALSE}
ncomps.csf_mxd35.plsda <- 3 # Number of component for PLS-DA 
set.seed(20220620) ## To replicate the results
start.time <- Sys.time() # Start the clock!
nsim <- 100 ## Number of repetitions

csf_mxd35_GPS.plsda.test.preds <- vector(mode="list", length=nsim)
csf_mxd35_GPS.plsda.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  train.idx <- partition(csf_mxd$Type, p = c(train=0.8, test=0.2)) ## library(SpliTools)
  
  train <- csf_mxd[train.idx$train,]
  test  <- csf_mxd[train.idx$test,]
  
  train.X <- as.matrix(train[,-1])
  train.Y<- as.numeric(as.character(train$Type))
  
  test.X <- as.matrix(test[,-1])
  test.Y <- as.numeric(as.character(test$Type))
  
  ## Select 35 Metabolites by GPS
  #m35_GPS <- gps(train.Y, train.X, family='binomial', maxvar=35, maxcmp=4)
  m35_GPS <- gps(train.Y, train.X, family='binomial', maxvar=35, maxcmp=5)
  
  ## PLS-DA with 35 Predictors Selected by GPS
  csf_mxd35_GPS.plsda.train <- plsda(train.X[,m35_GPS], train.Y, ncomp=ncomps.csf_mxd35.plsda)
  
  ## Predict
  csf_mxd35_GPS.plsda.test.preds[[j]] <- predict(csf_mxd35_GPS.plsda.train, test.X[,m35_GPS],
                                               dist="mahalanobis.dist")$predict[,,ncomps.csf_mxd35.plsda][,2]
  csf_mxd35_GPS.plsda.test.labels[[j]]<- test.Y
}

csf_mxd35_GPS.plsda.ROCR = ROCR::prediction(csf_mxd35_GPS.plsda.test.preds, csf_mxd35_GPS.plsda.test.labels)
csf_mxd35_GPS.plsda.roc = performance(csf_mxd35_GPS.plsda.ROCR,"tpr","fpr")
csf_mxd35_GPS.plsda.auc = performance(csf_mxd35_GPS.plsda.ROCR, measure="auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```




```{r}
## ROC curve and AUC
mar.default <- c(5,4,4,2) + 0.1
par(mar = mar.default + c(0, 2, 0, 0)) 
plot(csf_mxd35_GPS.plsda.roc,col="grey50", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_mxd35_GPS.plsda.roc, lwd=3, avg="vertical", cex.lab=1.5, cex.main=1.5, xlab="", ylab="",
     main="PLS-DA with 35 Predictors Selected by GPS")
text(x=0.8, y=0.1,paste("AUC=", round(median(unlist(csf_mxd35_GPS.plsda.auc@y.values)),2)),cex=1.5) # Show AUC
abline(a=0, b=1) 
par(new=FALSE)
```



```{r}
round(median(unlist(csf_mxd35_GPS.plsda.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_mxd35_GPS.plsda.auc@y.values)),2)/10  ## standard error
```



###   PLS-DA with 35 Predictors Selected by SIS

Tune the number of components for the PLS-DA

```{r}

set.seed(20220620) # for reproducibility, only when the `cpus' argument is not used
start.time <- Sys.time() # Start the clock!

train.idx <- partition(csf_mxd$Type, p = c(train=0.8, test=0.2)) ## library(SpliTools)
  
  train <- csf_mxd[train.idx$train,]
  test  <- csf_mxd[train.idx$test,]
  
  train.X <- as.matrix(train[,-1])
  train.Y<- as.numeric(as.character(train$Type))
  
  m35_SIS <- SIS(train.X, train.Y, family='binomial', iter=FALSE, tune='cv', nsis=35)$ix0
  m35_SIS_data=data.frame(Type=train$Type, train.X[, m35_SIS])
  
  m35_SIS.plsda <- plsda(m35_SIS_data[,-1], m35_SIS_data$Type, ncomp=10)  # set ncomp to 10


m35_SIS.plsda.perf <- perf(m35_SIS.plsda, validation="Mfold", folds=5,
                         progressBar=FALSE, auc=TRUE, nrepeat=50)
                         
ncomps.m35_SIS.plsda <- m35_SIS.plsda.perf$choice.ncomp
ncomps.m35_SIS.plsda

m35_SIS.plsda.perf$error.rate

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken) 
```


Run the model

```{r, message=FALSE, warning=FALSE}
set.seed(20220620) ## To replicate the results
ncomps.csf_mxd35.plsda <- 1  ## Number of component for PLS-DA 
start.time <- Sys.time() # Start the clock!
nsim <- 100 ## Number of repetitions

csf_mxd35_SIS.plsda.test.preds <- vector(mode="list", length=nsim)
csf_mxd35_SIS.plsda.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  train.idx <- partition(csf_mxd$Type, p=c(train=0.8, test=0.2)) ## library(SpliTools)
  
  train <- csf_mxd[train.idx$train,]
  test  <- csf_mxd[train.idx$test,]
  
  train.X = as.matrix(train[,-1])
  train.Y<- as.numeric(as.character(train$Type))
  
  test.X = as.matrix(test[,-1])
  test.Y = as.numeric(as.character(test$Type))
  
  # Use SIS to select variables
  #m35_SIS <- SIS(train.X, train.Y, family='binomial', iter=FALSE, tune='cv', nsis=35)$ix0
  m35_SIS <- SIS(train.X, train.Y, family='binomial', iter=FALSE, tune='cv', nsis=35)$ix0

  ## PLS-DA
  csf_mxd35_SIS.plsda.train <- plsda(train.X[,m35_SIS], train.Y, ncomp=ncomps.csf_mxd35.plsda)
  
  ## Predict
  csf_mxd35_SIS.plsda.test.preds[[j]] <- predict(csf_mxd35_SIS.plsda.train, test.X[,m35_SIS],
                                               dist="centroids.dist")$predict[,,ncomps.csf_mxd35.plsda][,2]
  csf_mxd35_SIS.plsda.test.labels[[j]]<- test.Y
}

csf_mxd35_SIS.plsda.ROCR = ROCR::prediction(csf_mxd35_SIS.plsda.test.preds, csf_mxd35_SIS.plsda.test.labels)
csf_mxd35_SIS.plsda.roc = performance(csf_mxd35_SIS.plsda.ROCR,"tpr","fpr")
csf_mxd35_SIS.plsda.auc = performance(csf_mxd35_SIS.plsda.ROCR, measure="auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```



```{r}
## ROC curve and AUC
mar.default <- c(5,4,4,2) + 0.1
par(mar = mar.default + c(0, 2, 0, 0)) 
plot(csf_mxd35_SIS.plsda.roc, col="grey50", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_mxd35_SIS.plsda.roc, lwd=3, avg="vertical", xlab="", ylab="",
      cex.lab=1.5, cex.main=1.5, main="PLS-DA with 35 Predictors Selected by SIS")
text(x=0.8, y=0.1,paste("AUC=", round(median(unlist(csf_mxd35_SIS.plsda.auc@y.values)),2)), cex=1.5)  # Show AUC
abline(a=0, b=1) 
par(new=FALSE)
```

```{r}
round(median(unlist(csf_mxd35_SIS.plsda.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_mxd35_SIS.plsda.auc@y.values)),2)/10  ## standard error
```



###  XGBoost-T with 35 Predictors Selected by GPS

** Tune the parameters for the XGBoost Tree **

```{r}
set.seed(20220620) # for reproducibility, only when the `cpus' argument is not used
start.time <- Sys.time() # Start the clock!

train.idx <- partition(csf_mxd$Type, p = c(train=0.8, test=0.2)) ## library(SpliTools)
  
train <- csf_mxd[train.idx$train,]
test  <- csf_mxd[train.idx$test,]
  
train.X <- as.matrix(train[,-1])
train.Y<- as.numeric(as.character(train$Type))
  
m35_GPS <- gps(train.Y, train.X, family='binomial', maxvar=35, maxcmp=5)

xgbtree_caret <- train(x=as.matrix(train.X[, m35_GPS]) ,
                   y=as.factor(train.Y),
                   method="xgbTree",
                   trControl=trainControl(method="repeatedcv",
                                          number=3,
                                          repeats=2,
                                          verboseIter=FALSE),
                   tuneGrid=expand.grid(nrounds=c(50,100,500),
                                        eta=c(0.01,0.05,0.1),
                                        max_depth=6,
                                        colsample_bytree=c(0.5, 1),
                                        subsample=c(0.5,1),
                                        gamma=0,
                                        min_child_weight=c(1,2)))
  
end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken) 

xgbtree_caret$bestTune ## Best parameters 
```



```{r}
set.seed(20220620) # for reproducibility here
start.time <- Sys.time() # Start the clock!
nsim <- 100    # number of simulations, could vary

csf_mxd35_GPS.xgbtree.test.preds <- vector(mode="list", length=nsim)
csf_mxd35_GPS.xgbtree.test.labels<- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  # Random Split into training and test data sets
  train.idx <- partition(csf_mxd$Type, p=c(train=0.8, test=0.2)) ## library(SpliTools)
  
  train <- csf_mxd[train.idx$train,]
  test <- csf_mxd[train.idx$test,]
  
  train.X <- as.matrix(train[,-1])
  train.Y <- as.numeric(as.character(train$Type))
  
  test.X = as.matrix(test[,-1])
  test.Y <- as.numeric(as.character(test$Type))
  
  ## Select 35 Metabolites by GPS
  m35_GPS <- gps(train.Y, train.X, family='binomial', maxvar=35, maxcmp=3)

  ## Convert the train and test data into xgb matrix type (required by XGboost package)
  xgb_train <- xgb.DMatrix(data=train.X[,m35_GPS], label=train.Y)
  xgb_test <- xgb.DMatrix(data=test.X[,m35_GPS], label=test.Y)
  
  ## Defining the parameters
  xgb_params <- list("objective"="binary:logistic",   ## logistic regression
                     "eval_metric"="logloss") ## multiclass logloss (used in classification)
  watchlist <- list(train=xgb_train, test=xgb_test)
  
  xgbc <- xgb.train(booster="gbtree",
                    params=xgb_params,
                    data=xgb_train,
                    max_depth=6,
                    eta=xgbtree_caret$bestTune$eta,
                    gamma=0,
                    nthread=1,
                    nrounds=xgbtree_caret$bestTune$nrounds, ## number of XGBoost iterations
                    watchlist=watchlist,
                    min_child_weight=xgbtree_caret$bestTune$min_child_weight, ##
                    subsample=xgbtree_caret$bestTune$subsample, 
                    colsample_bytree=xgbtree_caret$bestTune$colsample_bytree, ##
                    verbose=F) ## Print or not

  ### print(xgbc)
  
  ## The model is ready and we can predict the test data.
  ## Predicting test data (Hold-out)
  csf_mxd35_GPS.xgbtree.test.preds[[j]] <- as.numeric(predict(xgbc, xgb_test))
  csf_mxd35_GPS.xgbtree.test.labels[[j]] <- test.Y
}

csf_mxd35_GPS.xgbtree.ROCR <- ROCR::prediction(csf_mxd35_GPS.xgbtree.test.preds, csf_mxd35_GPS.xgbtree.test.labels)
csf_mxd35_GPS.xgbtree.roc <- performance(csf_mxd35_GPS.xgbtree.ROCR,"tpr","fpr")
csf_mxd35_GPS.xgbtree.auc <- performance(csf_mxd35_GPS.xgbtree.ROCR, measure = "auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```



```{r}
mar.default <- c(5,4,4,2) + 0.1
par(mar = mar.default + c(0, 2, 0, 0)) 
plot(csf_mxd35_GPS.xgbtree.roc, col="grey", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_mxd35_GPS.xgbtree.roc, lwd=3, avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="XGBoost-T with 35 Predictors Selected by GPS")
text(x=0.8, y=0.1,paste("AUC=", round(median(unlist(csf_mxd35_GPS.xgbtree.auc@y.values)),2)), cex=1.5) # AUC
abline(a=0, b=1) 
par(new=FALSE)
```

```{r}
round(median(unlist(csf_mxd35_GPS.xgbtree.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_mxd35_GPS.xgbtree.auc@y.values)),2)/10  ## standard error
```



###  XGBoost-T with 35 Predictors Selected by SIS

**Tune parameters **

```{r}
set.seed(20220620) # for reproducibility, only when the `cpus' argument is not used
start.time <- Sys.time() # Start the clock!

train.idx <- partition(csf_mxd$Type, p = c(train=0.8, test=0.2)) ## library(SpliTools)
  
train <- csf_mxd[train.idx$train,]
test  <- csf_mxd[train.idx$test,]
  
train.X <- as.matrix(train[,-1])
train.Y<- as.numeric(as.character(train$Type))
  
m35_SIS <- SIS(train.X, train.Y, family='binomial', iter=FALSE, tune='cv', nsis=35)$ix0
  
xgbtree_caret <- train(x=as.matrix(train.X[, m35_SIS]) ,
                   y=as.factor(train.Y),
                   method="xgbTree",
                   trControl=trainControl(method="repeatedcv",
                                          number=3,
                                          repeats=2,
                                          verboseIter=FALSE),
                   tuneGrid=expand.grid(nrounds=c(50,100,500),
                                        eta=c(0.01,0.05,0.1),
                                        max_depth=6,
                                        colsample_bytree=c(0.5,1),
                                        subsample=c(0.5,1),
                                        gamma=0,
                                        min_child_weight=c(1,2)))
  
end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken) 

xgbtree_caret$bestTune ## Best parameters 
```


```{r, message=FALSE, warning=FALSE}
set.seed(20220620) # for reproducibility here
start.time <- Sys.time() # Start the clock!
nsim <- 100    # number of simulations, could vary

csf_mxd35_SIS.xgbtree.test.preds <- vector(mode="list", length=nsim)
csf_mxd35_SIS.xgbtree.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  # Randomly split into training & test data sets
  train.idx <- partition(csf_mxd$Type, p=c(train=0.8, test=0.2)) ## library(SpliTools)

  train <- csf_mxd[train.idx$train,]
  test  <- csf_mxd[train.idx$test,]
  
  train.X <- as.matrix(train[,-1])
  train.Y <- as.numeric(as.character(train$Type))
  
  test.X = as.matrix(test[,-1])
  test.Y <- as.numeric(as.character(test$Type))
  
  # Use SIS to select variables
  m35_SIS <- SIS(train.X, train.Y, family='binomial', iter=FALSE, tune='cv', nsis=35)$ix0

  ## Convert the train and test data into xgb matrix type (required by XGboost package)
  xgb_train <- xgb.DMatrix(data=train.X[,m35_SIS], label=train.Y)
  xgb_test <- xgb.DMatrix(data=test.X[,m35_SIS], label=test.Y)
  
  ## Defining the parameters
  xgb_params <- list("objective"="binary:logistic",  ## logistic regression
                     "eval_metric" = "logloss") ## multiclass logloss (used in classification)
  watchlist <- list(train=xgb_train, test=xgb_test)

  xgbc <- xgb.train(booster="gbtree",
                    params=xgb_params,
                    data=xgb_train,
                    max_depth=6,
                    eta=xgbtree_caret$bestTune$eta,
                    gamma=0,
                    nthread=1,
                    nrounds=xgbtree_caret$bestTune$nrounds, ## number of XGBoost iterations
                    watchlist=watchlist,
                    min_child_weight=xgbtree_caret$bestTune$min_child_weight, ##
                    subsample=xgbtree_caret$bestTune$subsample, 
                    colsample_bytree=xgbtree_caret$bestTune$colsample_bytree, ##
                    verbose=F) ## Print or not
 
  ## print(xgbc)
  
  ## The model is ready and we can predict the test data.
  ## Predicting test data (Hold-out)
  csf_mxd35_SIS.xgbtree.test.preds[[j]] <- as.numeric(predict(xgbc, xgb_test))
  csf_mxd35_SIS.xgbtree.test.labels[[j]] <- test.Y
}

csf_mxd35_SIS.xgbtree.ROCR <- ROCR::prediction(csf_mxd35_SIS.xgbtree.test.preds, csf_mxd35_SIS.xgbtree.test.labels)
csf_mxd35_SIS.xgbtree.roc <- performance(csf_mxd35_SIS.xgbtree.ROCR,"tpr","fpr")
csf_mxd35_SIS.xgbtree.auc <- performance(csf_mxd35_SIS.xgbtree.ROCR, measure="auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```




```{r}
mar.default <- c(5,4,4,2) + 0.1
par(mar = mar.default + c(0, 2, 0, 0)) 
plot(csf_mxd35_SIS.xgbtree.roc, col="grey", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_mxd35_SIS.xgbtree.roc, lwd=3, avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="XGBoost-T with 35 Predictors Selected by SIS")
text(x=0.8, y=0.1,paste("AUC=", round(median(unlist(csf_mxd35_SIS.xgbtree.auc@y.values)),2)), cex=1.5) # AUC
abline(a=0, b=1) 
par(new=FALSE)
```


```{r}
round(median(unlist(csf_mxd35_SIS.xgbtree.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_mxd35_SIS.xgbtree.auc@y.values)),2)/10  ## standard error
```




### XGBoost-L with 35 Predictors Selected by GPS


**Tune parameters**

```{r, message=FALSE, warning=FALSE}
set.seed(20220620)
start.time <- Sys.time() # Start the clock!

train.idx <- partition(csf_mxd$Type, p = c(train=0.8, test=0.2)) ## library(SpliTools)
  
train <- csf_mxd[train.idx$train,]
test  <- csf_mxd[train.idx$test,]
  
train.X <- as.matrix(train[,-1])
train.Y<- as.numeric(as.character(train$Type))
  
m35_GPS <- gps(train.Y, train.X, family='binomial', maxvar=35, maxcmp=5)
xgblin_caret <- train(x=as.matrix(train.X) ,
                   y=as.factor(train.Y),
                   method="xgbLinear",
                   trControl=trainControl(method="repeatedcv",
                                          number=5,
                                          repeats=3,
                                          verboseIter=FALSE),
                   tuneGrid=expand.grid(nrounds=c(1,2,50),
                                        lambda =0,
                                        alpha = 1,
                                        eta=0.3
                                          ))

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```

The best parameters

```{r}
xgblin_caret$bestTune 

```


```{r}
xgblin_tuned <- xgb.train(data=xgb_train, booster="gblinear",
                        params=xgb_params,
                        nrounds=xgblin_caret$bestTune$nrounds,
                        verbose=0,
                        watchlist=watchlist)
```



```{r, message=FALSE, warning=FALSE}
set.seed(20220620) # for reproducibility here
start.time <- Sys.time() # Start the clock!
nsim <- 100    # number of simulations, could vary

csf_mxd35_GPS.xgblin.test.preds <- vector(mode="list", length=nsim)
csf_mxd35_GPS.xgblin.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  # Randomly split into training & test data sets
  train.idx <- partition(csf_mxd$Type, p=c(train=0.8, test=0.2)) ## library(SpliTools)
  
  train <- csf_mxd[train.idx$train,]
  test  <- csf_mxd[train.idx$test,]
  
  train.X <- as.matrix(train[,-1])
  train.Y <- as.numeric(as.character(train$Type))
  
  test.X <- as.matrix(test[,-1])
  test.Y <- as.numeric(as.character(test$Type))
  
  ## Select 35 Metabolites by GPS
  m35_GPS <- gps(train.Y, train.X, family='binomial', maxvar=35, maxcmp=3)

  ## Convert the train and test data into xgb matrix type (required by XGboost package)
  xgb_train <- xgb.DMatrix(data=train.X[,m35_GPS], label=train.Y)
  xgb_test <- xgb.DMatrix(data=test.X[,m35_GPS], label=test.Y)
  
  ## Defining the parameters
  xgb_params <- list("objective"="binary:logistic",  ## logistic regression
                     "eval_metric" = "logloss") ## multiclass logloss (used in classification)
  watchlist <- list(train=xgb_train, test=xgb_test)

  xgbc <- xgb.train(data=xgb_train, booster="gblinear",
                    params=xgb_params,
                    nrounds=xgblin_caret$bestTune$nrounds ,
                    verbose=0,
                    watchlist=watchlist)
  ### print(xgbc)
  
  ## The model is ready and we can predict the test data.
  ## Predicting test data (Hold-out)
  csf_mxd35_GPS.xgblin.test.preds[[j]] <- as.numeric(predict(xgbc, xgb_test))
  csf_mxd35_GPS.xgblin.test.labels[[j]] <- test.Y
  
}

csf_mxd35_GPS.xgblin.ROCR <- ROCR::prediction(csf_mxd35_GPS.xgblin.test.preds, csf_mxd35_GPS.xgblin.test.labels)
csf_mxd35_GPS.xgblin.roc <- performance(csf_mxd35_GPS.xgblin.ROCR,"tpr","fpr")
csf_mxd35_GPS.xgblin.auc <- performance(csf_mxd35_GPS.xgblin.ROCR, measure="auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```


```{r}
mar.default <- c(5,4,4,2) + 0.1
par(mar = mar.default + c(0, 2, 0, 0)) 
plot(csf_mxd35_GPS.xgblin.roc, col="grey", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_mxd35_GPS.xgblin.roc, lwd=3, avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="XGBoost-L with 35 Predictors Selected by GPS")
text(x=0.8, y=0.1, paste("AUC=", round(median(unlist(csf_mxd35_GPS.xgblin.auc@y.values)),2)), cex=1.5) # AUC
abline(a=0, b=1) 
par(new=FALSE)
```


```{r}
round(median(unlist(csf_mxd35_GPS.xgblin.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_mxd35_GPS.xgblin.auc@y.values)),2)/10  ## standard error
```




### XGBoost-L with 35 Predictors Selected by SIS

**Tune parameters**

```{r, message=FALSE, warning=FALSE}
set.seed(20220620)
start.time <- Sys.time() # Start the clock!

train.idx <- partition(csf_mxd$Type, p = c(train=0.8, test=0.2)) ## library(SpliTools)
  
train <- csf_mxd[train.idx$train,]
test  <- csf_mxd[train.idx$test,]
  
train.X <- as.matrix(train[,-1])
train.Y<- as.numeric(as.character(train$Type))
  
m35_SIS <- SIS(train.X, train.Y, family='binomial', iter=FALSE, tune='cv', nsis=35)$ix0

xgblin_caret <- train(x=as.matrix(train.X) ,
                   y=as.factor(train.Y),
                   method="xgbLinear",
                   trControl=trainControl(method="repeatedcv",
                                          number=5,
                                          repeats=10,
                                          verboseIter=FALSE),
                   tuneGrid=expand.grid(nrounds=c(1,2,5),
                                        lambda =0,
                                        alpha = 1,
                                        eta=0.3))

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```

The best parameters

```{r}
xgblin_caret$bestTune 
```


**Final model**

```{r, message=FALSE, warning=FALSE}
set.seed(20220620) # for reproducibility here
start.time <- Sys.time() # Start the clock!
nsim <- 100    # number of simulations, could vary

csf_mxd35_SIS.xgblin.test.preds <- vector(mode="list", length=nsim)
csf_mxd35_SIS.xgblin.test.labels <- vector(mode="list", length=nsim)

for(j in 1:nsim) {
  # Randomly split into training & test data sets
  train.idx <- partition(csf_mxd$Type, p=c(train=0.8, test=0.2)) ## library(SpliTools)
  
  train <- csf_mxd[train.idx$train,]
  test  <- csf_mxd[train.idx$test,]
  
  train.X <- as.matrix(train[,-1])
  train.Y<- as.numeric(as.character(train$Type))
  
  test.X <- as.matrix(test[,-1])
  test.Y <- as.numeric(as.character(test$Type))
  
  # Use SIS to select variables
  m35_SIS <- SIS(train.X, train.Y, family='binomial', iter=FALSE, tune='cv', nsis=35)$ix0

  ## Convert the train and test data into xgb matrix type (required by XGboost package)
  xgb_train <- xgb.DMatrix(data=train.X, label=train.Y)
  xgb_test <- xgb.DMatrix(data=test.X, label=test.Y)
  
  ## Defining the parameters
  xgb_params <- list("objective"="binary:logistic",  ## logistic regression
                     "eval_metric" = "logloss") ## multiclass logloss (used in classification)
  watchlist <- list(train=xgb_train, test=xgb_test)
  
  xgbc <- xgb.train(data=xgb_train, booster="gblinear",
                    params=xgb_params,
                    nrounds=xgblin_caret$bestTune$nrounds ,
                    verbose = 0,
                    watchlist=watchlist)
  ### print(xgbc)
  
  ## The model is ready and we can predict the test data.
  ## Predicting test data (Hold-out)
  csf_mxd35_SIS.xgblin.test.preds[[j]] <- as.numeric(predict(xgbc, xgb_test))
  csf_mxd35_SIS.xgblin.test.labels[[j]] <- test.Y
}

csf_mxd35_SIS.xgblin.ROCR = ROCR::prediction(csf_mxd35_SIS.xgblin.test.preds, csf_mxd35_SIS.xgblin.test.labels)
csf_mxd35_SIS.xgblin.roc = performance(csf_mxd35_SIS.xgblin.ROCR,"tpr","fpr")
csf_mxd35_SIS.xgblin.auc = performance(csf_mxd35_SIS.xgblin.ROCR, measure="auc")

end.time <- Sys.time() # Stop the clock!
time.taken <- end.time - start.time
print(time.taken)
```



```{r}
mar.default <- c(5,4,4,2) + 0.1
par(mar = mar.default + c(0, 2, 0, 0)) 
plot(csf_mxd35_SIS.xgblin.roc, col="grey", lty=3, xlab="False Positive Rate", ylab="True Positive Rate", cex.lab=1.5)
par(new=TRUE)
plot(csf_mxd35_SIS.xgblin.roc, lwd=3, avg="vertical", xlab="", ylab="",
     cex.lab=1.5, cex.main=1.5, main="XGBoost-L with 35 Predictors Selected by SIS")
text(x = 0.8, y = 0.1,paste("AUC=", round(median(unlist(csf_mxd35_SIS.xgblin.auc@y.values)),2)),cex = 1.5)  # Add text element
#mtext("B", side=3, line=1, cex=1.5, adj=-0.05)
abline(a=0, b=1) 
par(new=FALSE)
```



```{r}
round(median(unlist(csf_mxd35_SIS.xgblin.auc@y.values)),2) ## median AUC
round(sd(unlist(csf_mxd35_SIS.xgblin.auc@y.values)),2)/10  ## standard error
```



* AUC summary table for the **csf_mxd** study.

```{r}
Summary_AUC <- data.frame(Variable =c('All','35 with GPS','35 with SIS'),
               Logistic=c(0.56, '0.69','0.64'), PLSDA=c(0.71, 0.73, 0.72),
               XGBoostT=c(0.73, 0.73, 0.72), XGBoostL=c(0.70, 0.71, 0.69)
                                       )
Summary_AUC
```

